Global Navigation Satellite System Signal Decomposition and Parameterization Algorithm

ABSTRACT

A method and apparatus is provided for intra-PIT signal decomposition of a signal received with RF front end hardware. The method begins by aligning a signal received by RF front end hardware into integer multiples of a duration of a pseudorandom noise code sequence. A search grid is computed based on an integer multiple of the aligned signal. A plurality of initial ray parameters associated with the computed search grid is coarsely estimated. Using the coarsely estimated plurality of initial ray parameters, a fine estimation of the plurality of initial ray parameters is initiated utilizing stochastic search and optimization techniques. A stopping criteria statistic is computed by comparing a peak power of the search grid with a noise power present in the search grid. Finally, in response to determining the stopping criteria statistic being less than a stopping criteria threshold, processing a next integer multiple of the aligned signal.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of and priority to U.S. ProvisionalApplication Ser. No. 61/660,852, entitled “Global Navigation SatelliteSystem Signal Decomposition and Parameterization Algorithm,” filed onJun. 18, 2012, the entirety of which is incorporated by referenceherein.

RIGHTS OF THE GOVERNMENT

The invention described herein may be manufactured and used by or forthe Government of the United States for all governmental purposeswithout the payment of any royalty.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention generally relates to Global Navigation SatelliteSystem (GNSS) signals and, more particularly, to signal decompositionfor direct and indirect rays present in the GNSS data.

2. Description of the Related Art

Multipath is typically a major source of error in contemporary GNSSsystems. Ray models are generally used to describe GPS signals that arereceived in the presence of multipath propagation effects and whiteGaussian noise. Specular reflections in particular are typicallydescribed using a ray model. A Multipath Estimating Delay Lock Loop(MEDLL) receiver is a prominent example of a receiver design employedspecifically to mitigate the effects of multipath. The MEDLL receivermakes explicit use of multipath parameter estimation and tracking toreduce multipath effects on receiver outputs. MEDLL works bysimultaneously estimating parameters for both line-of-sight andmultipath signals. There are simultaneously operating tracking loopsemployed to track each of the signals present, both line-of-sight andmultipath.

The mathematics in MEDLL has been described to be akin to performingnonlinear curve fits. A set of reference correlation vectors aregenerated, each of which has an amplitude, phase, and propagation delaythat would be associated with a signal present in received data. Maximumlikelihood estimation is then used to obtain multipath signal parameterestimates that minimize the error between received signals and estimatedsignals generated by the receiver. The number of distinct signalsgenerated within the receiver can be varied, depending on how manymultipath signals are detected. MEDLL is perhaps the most prevalentmultipath mitigation technique that can be found in GNSS literature inwhich multipath estimates are explicitly obtained. However, adisadvantage of the conventional MEDLL receiver is a feedback loop,which may improperly bias future loop outputs based on past results thathave been erroneously obtained.

Commonly used multipath models do not adequately represent highlycomplex GNSS signal environments encountered in urban or indoorenvironments or under foliage. Because of this, a mitigation techniquehas been developed and verified in the lab, only to be deployed for useunsuccessfully because of the model's lack of adherence to thecharacteristics of the true environment. This problem may be alleviatedthrough the use of a recording and playback system. A system such asthis is used to record real-world GNSS signals, and then replay theserecorded signals. This is distinct from a simulator that uses models toconstruct signals that demonstrate propagation effects that would bepresent in the channel. The advantage of a recording and playback systemis that there is no uncertainty regarding the validity of the effectsobserved in the recorded signal. The disadvantage of this system is thatno models are obtained that can be used in mathematically describing thechannel.

Accordingly, there is a need in the art for a modeling system thatoffers both the advantages of a recording and playback system and asystem that uses an analytic model to generate data without the feedbackeffects seen in a MEDLL receiver.

SUMMARY OF THE INVENTION

Embodiments of the invention provide a method for intra-PIT signaldecomposition of a signal received with RF front end hardware. Themethod begins by aligning a signal received by RF front end hardwareinto integer multiples of a duration of a pseudorandom noise codesequence. A search grid is computed based on the aligned signal. Aplurality of initial ray parameters associated with the computed searchgrid is coarsely estimated. Using the coarsely estimated plurality ofinitial ray parameters, a fine estimation of the plurality of initialray parameters is initiated utilizing stochastic search and optimizationtechniques. An estimate error is computed for the fine estimation of theplurality of initial ray parameters. Finally, in response to determininga reduction in the computed estimate error, the fine estimation of theplurality of initial ray parameters is removed from the computed searchgrid.

In some embodiments of the invention, the plurality of initial rayparameters may include amplitude, carrier frequency, Doppler offset,propagation delay, carrier phase, and combinations thereof. In someembodiments, the received signal is a GPS L1 C/A-coded signal.

Other embodiments of the invention provide an apparatus having a RFfront end hardware configured to receive, downconvert, and sample asignal. The apparatus further includes a memory and a processor. Aprogram code is resident in the memory and configured to be executed bythe processor for intra-PIT signal decomposition of the signal receivedby the RF front end hardware. The program code is further configured toalign the signal into integer multiples of a duration of a pseudorandomnoise code sequence, compute a search grid based in the aligned signal,coarsely estimate a plurality of initial ray parameters associated withthe computed search grid, use the coarsely estimated plurality ofinitial ray parameters to initiate fine estimation of the plurality ofinitial ray parameters utilizing stochastic search and optimizationtechniques, compute an estimate error for the fine estimation of theplurality of initial ray parameters, and in response to determining areduction in the computed estimate error, remove the fine estimation ofthe plurality of initial ray parameters from the computed search grid.

Alternate embodiments of the invention provide a method for intra-PITsignal decomposition of a signal received with RF front end hardware.The method begins by aligning a signal received by RF front end hardwareinto integer multiples of a duration of a pseudorandom noise codesequence. A search grid is computed based on an integer multiple of thealigned signal. A plurality of initial ray parameters associated withthe computed search grid is coarsely estimated. Using the coarselyestimated plurality of initial ray parameters, a fine estimation of theplurality of initial ray parameters is initiated utilizing stochasticsearch and optimization techniques. A stopping criteria statistic iscomputed by comparing a peak power of the search grid with a noise powerpresent in the search grid. Finally, in response to determining thestopping criteria statistic is less than a stopping criteria threshold,processing a next integer multiple of the aligned signal.

Additional objects, advantages, and novel features of the invention willbe set forth in part in the description which follows, and in part willbecome apparent to those skilled in the art upon examination of thefollowing or may be leaned by practice of the invention. The objects andadvantages of the invention may be realized and attained by means of theinstrumentalities and combinations particularly pointed out in theappended claims.

BRIEF DESCRIPTION OF THE DRAWINGS

The accompanying drawings, which are incorporated in and constitute apart of this specification, illustrate embodiments of the invention and,together with a general description of the invention given above, andthe detailed description given below, serve to explain the invention.

FIG. 1 is a top-level flow chart illustrating a signal decompositionalgorithm consistent with embodiments of the invention;

FIG. 2 is a flow chart illustrating a signal parameterization algorithmconsistent with embodiments of the invention;

FIG. 3 is a flow chart illustrating a search space decomposition processconsistent with embodiments of the invention;

FIG. 4A is an illustration of a magnitude of points in a search grid forsimulated received data;

FIG. 4B is an illustration of a related phase of the points in FIG. 4A;

FIG. 5A is an illustration of a magnitude of points in a search grid fora trial estimate;

FIG. 5B is an illustration of a related phase of the points in FIG. 5A;

FIG. 6A is an illustration of a magnitude of points in a search grid fora difference between a received signal and a trail signal;

FIG. 6B is an illustration of a related phase of the points in FIG. 6A;

FIG. 7 is an illustration of a set of various trial parameters generatedusing a simulated annealing process consistent with embodiments of theinvention;

FIG. 8 is an illustration of errors associated with various trialparameters illustrated in FIG. 7;

FIG. 9 is an illustration of a trajectory of errors associated with anaddition of parameter estimates in an exemplary embodiment of theinvention;

FIG. 10 is a histogram illustrating normalized amplitude estimate errordistribution for an exemplary embodiment of the invention;

FIG. 11 is a histogram illustrating carrier frequency offset estimateerror distribution for the exemplary embodiment of FIG. 10;

FIG. 12 is a histogram illustrating propagation path length offset errordistribution for the exemplary embodiment of FIG. 10;

FIG. 13 is a histogram illustrating carrier phase error distribution forthe exemplary embodiment of FIG. 10;

FIG. 14 is a table of mean amplitude error (normalized) for an exemplaryembodiment of the invention;

FIG. 15 is a table of standard deviation of the amplitude error(normalized) of the table in FIG. 14;

FIG. 16 is a table of mean frequency offset error for the exemplaryembodiment of FIG. 14;

FIG. 17 is a table of standard deviation of the frequency offset errorof the table in FIG. 16;

FIG. 18 is a table of mean path length offset error for the exemplaryembodiment of FIG. 14;

FIG. 19 is a table of standard deviation of the path length offset errorof the table in FIG. 18;

FIG. 20 is a table of mean carrier phase error for the exemplaryembodiment of FIG. 14;

FIG. 21 is a table of standard deviation f the mean carrier phase errorof the table in FIG. 20;

FIG. 22 is a top-level flow chart illustrating an alternate a signaldecomposition algorithm consistent with embodiments of the invention;

FIGS. 23A and 23B is a flow chart illustrating a signal parameterizationalgorithm consistent with the embodiment of FIG. 22;

FIG. 24 is a flow chart illustrating a search space decompositionprocess consistent with embodiment FIG. 22; and

FIG. 25 is an exemplary hardware and software environment consistentwith embodiments of the invention.

It should be understood that the appended drawings are not necessarilyto scale, presenting a somewhat simplified representation of variousfeatures illustrative of the basic principles of the invention. Thespecific design features of the sequence of operations as disclosedherein, including, for example, specific dimensions, orientations,locations, and shapes of various illustrated components, will bedetermined in part by the particular intended application and useenvironment. Certain features of the illustrated embodiments have beenenlarged or distorted relative to others to facilitate visualization andclear understanding. In particular, thin features may be thickened, forexample, for clarity or illustration.

DETAILED DESCRIPTION OF THE INVENTION

Multipath Estimating Delay Lock Loop (MEDLL) receivers are prominentexamples of a receiver design employed specifically to mitigate theeffects of multipath. The MEDLL receiver makes explicit use of multipathparameter estimation and tracking to reduce multipath effects onreceiver outputs. MEDLL works by simultaneously estimating parametersfor both line-of-sight and multipath signals. There are simultaneouslyoperating tracking loops employed to track each of the signals present,both line-of-sight and multipath.

Like the MEDLL receiver, the embodiments of the invention set out inmore detail below specifically obtain estimates of multipath, as well asline of sight (if available) signals. However, one advantage of thealgorithm in the embodiments of the invention over the conventionalMEDLL receiver is the absence of any feedback loop which may improperlybias future loop outputs based on past results that have beenerroneously obtained. This is because a search space decompositionalgorithm is executed independently on search space outputs computedfrom individual integration periods, similarly to the computations thatwould be performed in block processing of GNSS signals. The algorithm inthe embodiments of the invention is also able to properly characterizemore complex signal environments than MEDLL.

In a harsh signal environment where RF propagators, such as reflectorsand scatterers, are present in a relative proximity to the GNSS receiverantenna, interference that is generated by these propagators willdegrade the integrity of the line-of-sight signal unless the effect ofthe interference is mitigated. However, because there most likely is adifference in the time of arrival of interfering signals, relative tothe line of sight signal, the constituent multipath elements of thereceived signal can be distinctly observed, using the frequency/codephase search space. Since these signal elements can be separatelyobserved, these elements can then be individually parameterized using aglobal stochastic search and optimization method, such as simulatedannealing or potentially a genetic algorithm.

Turning now to the drawings, FIG. 1 illustrates a top-level flowchart 10of an intra-PIT (pre-detection integration time) signal decompositionalgorithm. This figure corresponds directly with a process flowchart 12illustrated in FIG. 2. The algorithm begins at block 14 of flowchart 10in FIG. 1. RF front end functions (such as frequency downconversion andsignal sampling are performed in block 16. Then, the signal is alignedin block 18 so that data may be processed by a decomposition algorithmin integer multiples of a duration of the pseudorandom noise (PRN) codesequence present in GPS L1 C/A-coded signals, i.e., integer multiples of1,000 msec. The decomposition algorithm, in some embodiments, may onlyprocess whole periods of the PRN code sequence. For example, in aparticular embodiment, a PIT-sized data vector is processed, where thefirst sample in the data vector is from the first chip in the PRNsequence, and the last sample in the vector is from the last chip in thePRN sequence. Other embodiments may utilize decomposition algorithmsthat process PRN code sequences using other data sizes. Once alignmenttakes place, received data that is sized to be of PIT-sized duration isinput for search grid computation and initial ray parameter estimationin block 20. The amplitude, the carrier frequency Doppler offset, thepropagation delay, and the carrier phase are coarsely estimated. Thesecoarse parameter estimates may then be used to initiate fine estimationof these four parameters for a ray using a stochastic search andoptimization technique (simulated annealing, for example) in block 22.The error of the parameter estimates is then computed in block 24, andthe portion of the search grid associated with the estimated parametersis removed from the search grid in block 26 (“Yes” branch of decisionblock 28). The remainder of the search grid is then input into theprocess to obtain another ray estimate at block 20. This continues untilthere are no more rays to estimate (“No” branch of decision block 28).The next PIT-sized data vector is then input for processing at block 30.This loop takes place until all data is processed and parameterized.

Referring again to FIG. 2, flowchart 12 outlines a process used toestimate the multipath-propagated rays received with a software receiverfor some embodiments of the invention. This flowchart 12 involves theexecution of several processes, each of which are set out in furtherdetail below. Generally, these processes are as follows: downconversionof received data to an intermediate frequency, sampling of data,buffering of code-aligned data, computation of the search grid fromreceived data, estimation of initial parameters, computation of theinitial amplitude estimate, decomposition of the search space,computation of the estimate search space, comparison of received andestimate search spaces, computation of the estimate error, output ofparameters to the ray database, determination of alternate search gridcomparison regional maxima to consider, and initial estimation ofalternate peak parameters.

For the description of the algorithm implemented in embodiments of theinvention, the following model will be used for the continuoustime-domain received signal r(t):

$\begin{matrix}{{r(t)} = {{\sum\limits_{m = 0}^{M{(t)}}\; {{A_{m}(t)}{D\left( {t - \tau_{m)}} \right)}{x\left( {t - \tau_{m)}} \right)}{\cos \left( {{2\; {\pi \left( {f_{L\; 1} + {f_{D_{m}}(t)}} \right)}t} + \varphi_{r_{m}}} \right)}}} + {\eta (t)}}} & (1)\end{matrix}$

where m is an index of a multipath ray and an LOS ray is a 0th multipathray (m=0 for the ray assumed to be LOS), M(t) is a number of multipathrays at time t, A_(m)(t) is a peak amplitude of a signal from the mthray at time t, D(t−τ_(m)) is a value of the navigation data messagesymbol at time t−τ_(m), τ_(m) is a propagation delay of the mth raybetween a transmitting SV and a receiver, x(t−τ_(m)) is a value of aC/A-coded PRN sequence chip value at time t−τ_(m), f_(L1) is a centerfrequency of a GPS L1 band, f_(D) _(m) (t) is a Doppler frequency offsetfrom the L1 center frequency at time t, φ_(r) _(m) is a carrier phase ofthe mth multipath ray, and η(t) is white Gaussian noise that is added tothe received LOS and multipath signals at time t. This model will beused to fundamentally describe the signal at entry into the receiverthroughout the remainder of this description. It should be understoodhowever, that this model is merely a design choice, chosen forillustration purposes, and that other models may be used in otherembodiments.

Turning now to flowchart 12 in FIG. 2, the process begins at block 32.Received data is down converted at block 34. This process may beexecuted in RF front end hardware of a software receiver in someembodiments, though other configurations of the RF front end are alsocontemplated. The model that is used to describe the output of thedownconversion process is presented in further detail as follows. Thisprocess accepts r(t) as an input and outputs a downconverted receivedsignal r_(IF)(t) whose band has a center frequency of f_(IF). There arefrequencies other than f_(L1) that are available for various GPS users,but for the purposes of this description of the multipath estimationalgorithm, only L1 will be considered. The equation used to describe theoutput of the downconversion process is as follows:

$\begin{matrix}{{r_{IF}(t)} = {{\sum\limits_{m = 0}^{M{(t)}}\; {{A_{m}(t)}{D\left( {t - \tau_{m)}} \right)}{x\left( {t - \tau_{m)}} \right)}{\cos \left( {{2\; {\pi \left( {f_{IF} + {f_{D_{m}}(t)}} \right)}t} + \varphi_{m}} \right)}}} + {\eta (t)}}} & (2)\end{matrix}$

where f_(IF) is the intermediate frequency of the receiver, and φ_(m) isthe carrier phase of the mth multipath ray at the input into the ADC. Itis assumed that r_(IF)(t) is filtered using an ideal bandpass filtercentered at f_(IF). Both SiGe and TRIGR receivers use an intermediatefrequency when outputting downconverted and sampled signals forprocessing by a user. This process within a software receiver istransparent to the user. Note at this point that the carrier phase ofthe multipath ray input to the ADC may be completely different from thecarrier phase of the original received signal, thus the distinctionbetween φ_(r) _(m) in (1) and φ_(m) in (2). This does not impact how thesignal is processed. Again, the equation above is a model that makes useof the assumption of ideal bandpass filtering in its development, and aswith the model for the received signal, other models describing theoutput of the downconversion process may be used in other embodiments ofthe invention.

This equation for r(t) sufficiently serves as the model for the TRIGRreceiver, but does not describe the model to be used for the SiGereceiver. This is because the TRIGR receiver only processes areal-valued received signal. The SiGe receiver, on the other hand,processes a simultaneous in-phase signal i(t) and quadrature signalq(t). The model to describe continuous time-domain data where in-phaseand quadrature signals are both present is as follows:

r _(IF)(t)=i(t)+jq(t)+η_(i)(t)+jη _(q)(t)  (3)

where j=√{square root over (−1)}. Note the presence of complex-valuednoise in this model. This is done to account for the presence ofzero-mean white Gaussian noise that is added independently to both thein-phase and quadrature signals. The noise on the in-phase signal andthe noise on the quadrature signal is independent of each other, andshare the same variance. The continuous time-domain in-phase andquadrature signals i(t) and q(t) may be modeled as follows:

$\begin{matrix}{{{i(t)} = {\sum\limits_{m = 0}^{M{(t)}}\; {{A_{m}(t)}{D\left( {t - \tau_{m)}} \right)}{x\left( {t - \tau_{m)}} \right)}{\cos \left( {{2{\pi \left( {f_{IF} + {f_{D_{m}}(t)}} \right)}t} + \varphi_{r_{m}}} \right)}}}}\mspace{79mu} {and}} & (4) \\{{q(t)} = {\sum\limits_{m = 0}^{M{(t)}}\; {{A_{m}(t)}{D\left( {t - \tau_{m)}} \right)}{x\left( {t - \tau_{m)}} \right)}{\sin \left( {{2{\pi \left( {f_{IF} + {f_{D_{m}}(t)}} \right)}t} + \varphi_{r_{m}}} \right)}}}} & (5)\end{matrix}$

r_(if)(t)ε

is now restated using a complex exponential instead of sine and cosinefunctions:

$\begin{matrix}{{r_{IF}(t)} = {{\sum\limits_{m = 0}^{M{(t)}}\; {{A_{m}(t)}{D\left( {t - \tau_{m)}} \right)}{x\left( {t - \tau_{m)}} \right)}{\exp \left( {j\left( {{2{\pi \left( {f_{IF} + {f_{D_{m}}(t)}} \right)}t} + \varphi_{m}} \right)} \right)}}} + {\eta_{iq}(t)}}} & (6) \\{\mspace{76mu} {where}\mspace{20mu}} & \; \\{\mspace{79mu} {{\eta_{iq}(t)} = {{\eta_{i}(t)} + {{j\eta}_{q}(t)}}}} & (7)\end{matrix}$

r_(if)(t)ε

is further restated by dividing the complex exponential term into twoseparate terms for the carrier frequency and the carrier phase asfollows:

$\begin{matrix}{{r_{IF}(t)} = {{\sum\limits_{m = 0}^{M{(t)}}\; {{A_{m}(t)}{D\left( {t - \tau_{m)}} \right)}{x\left( {t - \tau_{m)}} \right)}{\exp \left( {j\left( {2{\pi \left( {f_{IF} + {f_{D_{m}}(t)}} \right)}t} \right)} \right)}}} + {\exp \left( {j\varphi}_{r_{m}} \right)} + {\eta_{iq}(t)}}} & (8)\end{matrix}$

After the downconversion process, data is sampled at block 36. Thisprocess may also be executed in the front end hardware of the softwarereceiver. The model that is used to describe the output of the samplingprocess is set out in more detail as follows. This process acceptsr_(IF)(t) as an input and outputs a discrete time-domain received signalr[n]. For bookkeeping purposes, the PIT index value will be initializedby asserting that p=0. The equation used to express the sampling ofr_(IF)(t), which yields r[n], is as follows.

r[n]=r _(IF)(nT _(s)),nε

  (9)

where n is the sample index and T_(s) is the sampling period. Thisprocess within a software receiver may again be transparent to the user.This yields the following result:

$\begin{matrix}{{r\lbrack n\rbrack} = {\left( {\sum\limits_{m = 0}^{M{\lbrack n\rbrack}}\; {{A_{m}\lbrack n\rbrack}{D_{m}\lbrack n\rbrack}{x_{m}\lbrack n\rbrack}{f_{m}\lbrack n\rbrack}{\theta_{m}\lbrack n\rbrack}}} \right) + {\eta \lbrack n\rbrack}}} & (10) \\{{M\lbrack n\rbrack} = {M\left( {nT}_{s} \right)}} & (11) \\{{A_{m}\lbrack n\rbrack} = {A_{m}\left( {nT}_{s} \right)}} & (12) \\{{D_{m}\lbrack n\rbrack} = {D\left( {{nT}_{s} - \tau_{m}} \right)}} & (13) \\{{x_{m}\lbrack n\rbrack} = {x\left( {{nT}_{s} - \tau_{m}} \right)}} & (14) \\{{f_{m}\lbrack n\rbrack} = {\exp \left( {j\left( {2{\pi \left( {f_{IF} + {f_{D_{m}}\left( {nT}_{s} \right)}} \right)}{nT}_{s}} \right)} \right)}} & (15) \\{{\theta_{m}\lbrack n\rbrack} = {\exp \left( {j\varphi}_{m} \right)}} & (16) \\{{\eta \lbrack n\rbrack} = {\eta_{iq}\left( {nT}_{s} \right)}} & (17)\end{matrix}$

Note the use of notation where only the discrete-time index is used (asin M[n]), rather than the use of the sampling period multiplied by thediscrete-time index (as in M(nT_(s))). This is to make use of thenotation convention to indicate that these equations, as well assubsequent equations that make use of the discrete-time index alone inbrackets, compose vectors of values over discrete time instances.

After sampling, the Code-aligned data is buffered in block 38. Thisprocess accepts r[n] and p as inputs and outputs a buffered vector ofreceived data r_(p) of size N_(p). The code-aligning buffer processesthe discrete time-domain received data r[n], where r[n]ε{

,

}, and divides r[n] into vectors of size N_(p). Therefore, r_(p) isexpressed as follows:

r _(p) =[r[n]] _(n=pN) _(p) ^(N) ^(p) ^(−1+pN) ^(p) ,pε

  (18)

However, this can be restated as:

$\begin{matrix}{r_{p} = \left\lbrack {{\sum\limits_{m = 0}^{M{\lbrack n\rbrack}}\; {{A_{m}\lbrack n\rbrack}{D_{m}\lbrack n\rbrack}{x_{m}\lbrack n\rbrack}{f_{m}\lbrack n\rbrack}{\theta_{m}\lbrack n\rbrack}}} + {\eta \lbrack n\rbrack}} \right\rbrack_{n = {pN}_{p}}^{N_{p} - 1 + {pN}_{p}}} & (19)\end{matrix}$

Because of this, (18) can be restated as a sum of individual vectors asfollows:

$\begin{matrix}{{r_{p} = {{\sum\limits_{m = 0}^{M_{p}}\; r_{p_{m}}} + \eta_{p}}}{where}} & (20) \\{{r_{p_{m}} = \left\lbrack {{A_{m}\lbrack n\rbrack}{D_{m}\lbrack n\rbrack}{x_{m}\lbrack n\rbrack}{f_{m}\lbrack n\rbrack}{\theta_{m}\lbrack n\rbrack}} \right\rbrack_{n = {pN}_{p}}^{N_{p} - 1 + {pN}_{p}}}{and}} & (21) \\{\eta_{p} = \left\lbrack {\eta \lbrack n\rbrack} \right\rbrack_{n = {pN}_{p}}^{N_{p} - 1 + {pN}_{p}}} & (22)\end{matrix}$

Note this expression is constrained by making the assumption that thenumber of rays is fixed within one PIT, thusM_(p)=M[n]∀nε{pN_(p),pN_(p)+1, . . . , N_(p−1)+pN_(p)}. (21) may berestated as follows:

r _(p) _(m) =A _(p) _(m) D _(p) _(m) x _(p) _(m) ⊙f _(p) _(m) θ_(p) _(m)  (23)

where ⊙ denotes a Hadamard product (element-by-element arraymultiplication). The amplitude A_(p) _(m) is constrained by assumingthat the amplitude for each ray is fixed within one PIT. The sameconstraining assumption is made for a navigation data message bit valueD_(p) _(m) and the complex exponential carrier phase term θ_(p) _(m) .

A _(p) _(m) =A _(m) [n]∀nε{pN _(p) ,pN _(p)+1, . . . ,N _(p−1) +pN_(p)}  (24)

D _(p) _(m) =D _(m) [n]∀nε{pN _(p) ,pN _(p)+1, . . . ,N _(p−1) +pN_(p)}  (25)

x _(p) =[x[n]] _(n=pN) _(p) ^(N) ^(p) ^(−1+pN) ^(p)   (26)

f _(p) =[f[n]] _(n=pN) _(p) ^(N) ^(p) ^(−1+pN) ^(p)   (27)

θ_(p) _(m) =θ_(m) [n]∀nε{pN _(p) ,pN _(p)+1, . . . ,N _(p−1) +pN_(p)}  (28)

Though not expressed in the equations above, it is assumed that theDoppler frequency offset f_(D) _(m) for each ray, as used in (15), isfixed within one PIT.

Given the C/A-coded PRN sequence vector b (of size N_(c)=1023) that ismodulated within a GPS waveform to spread the signal spectrum, alignmentof the C/A-coded PRN sequence from one PIT p to the next will be forced.This may be accomplished by simply ensuring that the most prominent rayreceived in r_(p) _(m) , which is assumed to be the LOS ray (denotedwith m=0), is output from the code-aligning buffer with the followingconstraint in place:

x _(p) ₀ =[b ₀ , . . . , b _(N) _(c) ⁻¹ ]∀p  (29)

This means that the first element of the coding vector x_(p) _(m) forthe LOS ray would contain the value of the first chip in the C/A-codedPRN sequence vector, and the last element in x_(p) _(m) would containthe value of the last chip in the sequence vector. This may be difficultin practical processing of data, but it will be assumed that this is thecase for this illustrated embodiment.

This aligning function may be performed by simply processing thereceived GPS data using a simple software receiver, and then trackingpoints in the signal where x[n]=b₀, while x[n−1]=b_(N) _(c) ⁻¹.

This function may be performed the same way regardless of the value ofPIT p.

After alignment, computation of a search grid from received data may beperformed in block 40. This process accepts r_(p) and N_(p) as inputsand outputs a search space R_(p) for the received signal vector r_(p).The search space may also be referred to here as a search grid. Becausethe same algorithm will be used for computing the search grid in severalinstances, with variation potentially in the propagation delay (oralternately path length offset) window size τ_(window) (or alternatelyδd_(window)), the frequency window size f_(window), the propagationdelay resolution (or alternately path length offset resolution) Δτ (oralternately Δ(δd)) and the frequency resolution Δf, the “grid” operatorwill be used to describe the search space operation, which yields anarbitrary search grid S:

S=grid(s,τ _(window) ,f _(window),Δτ_(input) ,Δf _(input))  (30)

where s is a vector of PIT size that is input into the search spacecomputation operation. The input value Δf_(input) is a desired frequencyresolution, as defined by a user. The input value Δτ_(input) is adesired propagation delay resolution, as defined by a user. The actualresolution of the search space, in both propagation delay and frequency,will be at least as small as the values input by the user forresolution.

Before detailing how the search grid is computed for the processoutlined, it is asserted that the search grid R_(p) for received PITsized vector r_(p) is computed as follows:

R _(p)=grid(r _(p),τ_(window) ,f _(window),Δτ_(input) ,Δf_(input))  (31)

where τ_(window)=2T_(c),f_(window)=10,000 kHz, Δτ_(input)=T_(s), andΔf_(input)=10 Hz. The chip period T_(c) is equal to the duration of onechip in the PRN sequence. The sample period T_(s)=1/f_(s), where f_(s)is the sample frequency.

The propagation delay window size τ_(window) is a user-dictated designchoice. The choice of τ_(window)=2T_(c) is made for this process inorder to obtain estimates of the multipath rays adjacent to the LOS ray,in order to characterize the channel for the local environment aroundthe receiver. Multipath much more prominently affects receiverprocessing when a multipath ray is received with a propagation delay ofless than T_(c), relative to the LOS ray. Since the delay between thepeak of a ray in the search space and the point where the ray isnegligible in power is T_(c), the sum of T_(c)+T_(c) is used for thewindow size. However, this process can easily be executed for anyarbitrary choice of propagation delay without requiring any changes tothe algorithm. The conversion between τ_(window) and the path lengthoffset window size δd_(window) is as follows:

δd _(window) =cτ _(window)  (32)

where c is the speed of light.

The frequency window size f_(window) is also a user-dictated designchoice. The choice of f_(window)=10,000 kHz is made for this illustratedembodiment because, as can be found in contemporary literature, Dopplershifts of up to f_(window)=10,000 kHz can be observed in receivedsignals, though other window sizes may also be used with otherembodiments of the invention.

The propagation delay resolution Δτ_(input) is a user-dictated designchoice. The choice of Δτ_(input)=T_(s) is made because this value of Δτcorresponds with the actual resolution afforded by the sample rate ofthe receiver. Other propagation rates may also be used to accommodatethe resolution of other receivers. The conversion betweenΔτ_(input)=T_(s) and the path length offset resolution Δ(δd) is asfollows:

Δ(δd)=cΔτ  (33)

The propagation delay resolution Δτ that is actually used in computationof the search grid is to always be at least as small as the resolutionspecified by the user. The equation for calculation of Δτ is as follows:

$\begin{matrix}{{{\Delta\tau} = \frac{2\tau_{window}}{N_{\tau} - 1}}{where}} & (34) \\{N_{\tau} = {{2\left\lceil {\max \left( {{f_{s}\tau_{window}},\frac{\tau_{window}}{{\Delta\tau}_{input}}} \right)} \right\rceil} + 1}} & (35)\end{matrix}$

The frequency resolution Δf that is actually used in computation of thesearch grid is to always exceed the resolution dictated by the user. Thesmaller the resolution, the better the initial estimate of the Dopplerfrequency offset that can be made for each ray in the PIT, thus theresolution choice of Δf_(input)=10 Hz. The search grid may be computedusing several methods, generally making use of a fast Fourier transform(FFT). Therefore, the potential use of the FFT will be considered whencomputing Δf, but will not be considered when presenting the generalform for the search grid in (41) below. The equation for calculation ofΔf is as follows:

$\begin{matrix}{{{\Delta \; f} = \frac{f_{s}}{N_{fft}}}{where}} & (36) \\{N_{fft} = {\max\left( {N_{s},2^{({\lceil{{\log \;}_{2}{(\frac{f_{s}}{\Delta \; f_{input}})}}\rceil})}} \right)}} & (37)\end{matrix}$

where N_(s) is the size of the data vectors input to the “grid”operation and ┌┐ indicates a ceiling rounding operation (rounding up tothe nearest integer on the real number line more than the valuecontained in the modified brackets). Note that the function expressedwith this equation is designed to make use of the computationalefficiency afforded the FFT function when input vectors to an FFTfunction are of a length equal to a power of two.

The fast Fourier transform algorithm “fft” is known in the art and willnot be further described, other than to say that the following operatoris used to denote the FFT operation, where s_(f) is the vector output ofthe FFT operation:

s _(f) =fft(s,N _(fft))_(fε[f) _(low) _(,f) _(high) _(])  (38)

where N_(fft) is the number of points in the input vector in the FFT,and is as expressed in (37). The expression fε[f_(low), f_(high)]denotes the range of frequencies, between f_(Low) and f_(high), forwhich data will be output in s_(f). Put another way, the vector outputfrom an FFT operation is associated with a range of frequencies asfollows: fε[−f_(s)/2, f_(s)/2]. The subscript to the fft operator isused to denote that the range of frequency-domain points output in s_(f)is limited to [f_(low),f_(high)], rather than the entire range thatwould be output were the FFT not constrained. The number of points inthe FFT whose frequencies are in the range fε[f_(low),f_(high)] isdenoted N_(f). The following vector γ is defined to be the frequencyvalues for each column in the search space:

γ=[γ_(z)]_(z=0) ^(N) ^(f) ⁻¹,γ_(z) ε[f _(low) ,f _(high)]  (39)

The code delay vector x_(s) for code delay τ_(s) is defined as follows:

$\begin{matrix}{{{x_{s}\left( \tau_{s} \right)} = \left\lbrack {b\left\lbrack \left\lfloor {N_{c}\left( {\frac{{\omega \; T_{s}} + \tau_{s}}{N_{c}T_{c}} - \left\lfloor \frac{{\omega \; T_{s}} + \tau_{s}}{N_{c}T_{c}} \right\rfloor} \right)} \right\rfloor \right\rbrack} \right\rbrack_{\omega = 0}^{N_{s} - 1}},{b\mspace{11mu} \varepsilon \left\{ {{- 1},{+ 1}} \right\}}} & (40)\end{matrix}$

where τ_(s) is an arbitrary propagation delay input into the expression,b[] is an element of the PRN sequence vector b designated by thenatural number inside the brackets, └┘ indicates a floor roundingoperation (rounding down to the nearest integer on the real number lineless than the value contained in the modified brackets), N_(c) is thenumber of chips in the PRN sequence, T_(c) is the chip duration (soT_(c)=0.001/N_(c)=0.001/1023 seconds for GPS C/A code), and N_(s) is asused in (37).

The search space computation operator “grid” and the output of the“grid” operator, S, is defined as follows, where the “grid” operator forS is first presented in (30) to generally relate the “grid” operator to(31):

$\begin{matrix}{S = {{{grid}\left( {s,\tau_{window},f_{window},{\Delta \; \tau_{input}},{\Delta \; f_{input}}} \right)} = \begin{bmatrix}\begin{matrix}{{S\left( {\beta_{0},\gamma_{0}} \right)}\;} \\{{S\left( {\beta_{1},\gamma_{0}} \right)}\;}\end{matrix} & \begin{matrix}{S\left( {\beta_{0},\gamma_{1}} \right)} \\{S\left( {\beta_{1},\gamma_{1}} \right)}\end{matrix} & \cdots & \begin{matrix}{S\left( {\beta_{0},\gamma_{N_{\tau} - 1}} \right)} \\{S\left( {\beta_{1},\gamma_{N_{\tau} - 1}} \right)}\end{matrix} \\\vdots & \; & \ddots & \vdots \\{{S\left( {\beta_{N_{\tau} - 1},\gamma_{0}} \right)}\;} & {S\left( {\beta_{N_{\tau} - 1},\gamma_{1}} \right)} & \cdots & {S\left( {\beta_{N_{\tau} - 1},\gamma_{N_{\tau} - 1}} \right)}\end{bmatrix}_{N_{\tau} \times N_{f}}}} & (41)\end{matrix}$

where Nτ is the size of the vector of code delays considered in S, thusthe number of rows in S, and is equal to the size of the followingvector β defining the code delay values for each row in the searchspace:

$\begin{matrix}{\beta = \left\lbrack {{- \tau_{window}},{{- \tau_{window}} + {\Delta\tau}},{{- \tau_{window}} + {2\; {\Delta\tau}}},\ldots \mspace{14mu},\tau_{window}} \right\rbrack^{T}} & (42) \\{\mspace{14mu} {= \left\lbrack \left\lbrack \beta_{y} \right\rbrack_{y = 0}^{N_{\tau} - 1} \right\rbrack^{T}}} & (43)\end{matrix}$

To compute S(β_(y),γ_(z)), the following equation is used:

$\begin{matrix}{\mspace{79mu} {{{S\left( {\beta_{y},\gamma_{z}} \right)} = {{S_{i}\left( {\beta_{y},\gamma_{z}} \right)} + {{jS}_{q}\left( {\beta_{y},\gamma_{z}} \right)}}}\mspace{79mu} {where}}} & (44) \\{{S_{i}\left( {\beta_{y},\gamma_{z}} \right)} = {\sum\limits_{\omega = 0}^{N_{s} - 1}\; \left( {\left( {{{{Re}\left( s_{\omega} \right)}{\cos \left( {\varphi_{s_{\omega}}\left( \gamma_{y} \right)} \right)}} + {{Im}\; \left( s_{\omega} \right){\sin \left( {\varphi_{s_{\omega}}\left( \gamma_{y} \right)} \right)}}} \right){x_{s_{\omega}}\left( \beta_{y} \right)}} \right)}} & (45) \\{{S_{q}\left( {\beta_{y},\gamma_{z}} \right)} = {\sum\limits_{\omega = 0}^{N_{s} - 1}\; \left( {\left( {{{{Im}\left( s_{\omega} \right)}{\cos \left( {\varphi_{s_{\omega}}\left( \gamma_{y} \right)} \right)}} - {R\; {e\left( s_{\omega} \right)}{\sin \left( {\varphi_{s_{\omega}}\left( \gamma_{y} \right)} \right)}}} \right){x_{s_{\omega}}\left( \beta_{y} \right)}} \right)}} & (46) \\{\mspace{79mu} {and}} & \; \\{\mspace{79mu} {{\varphi_{s_{\omega}}\left( \gamma_{y} \right)} = {2\; {\pi \left( {f_{IF} + \gamma_{y}} \right)}T_{s}\omega}}} & (47)\end{matrix}$

Returning to flowchart 12 in FIG. 2, estimation of initial parameters isthen performed in block 42. This process accepts R_(p) as an input andoutputs an initial estimate of the Doppler frequency offset {tilde over(f)}_(p) _(m) , an initial estimate of the propagation delay {tilde over(τ)}_(p) _(m) , and an initial estimate of the carrier phase {tilde over(φ)}_(p) _(m) . Additionally, in the case where the ray index m>0,{circumflex over (R)}_(p) _(m−1) is additionally input. {circumflex over(R)}_(p) _(m) is a search grid computed from the summation of discretetime-domain vectors associated with all previously obtained rayestimates. {circumflex over (R)}_(p) _(m) contains aggregated searchgrid information for m previously obtained ray estimates.

In the discussion of the computation of the search grid above, thenumber of rows in the search space matrix R_(p) is established to beN_(τ) and the number of columns in R_(p) to be N_(f). {circumflex over(R)}_(p) _(m−1) is additionally asserted in this process is of the samesize, and that the frequency and code phase values of element positionsshared between the two matrices are the same. The comparison matrixχ_(p) _(m) is now defined as follows:

$\begin{matrix}{\chi_{p_{m}} = \begin{Bmatrix}{{R_{p} - {\hat{R}}_{p_{m}}},} & {m > 0} \\{R_{p},} & {m = 0}\end{Bmatrix}} & (48)\end{matrix}$

Defining χ_(p) _(m) as shown in (48) simply provides a means by whichthe received search space is compared with the most recently computedsearch space of the sum of estimated multipath rays. As the number ofrays in {circumflex over (R)}_(p) _(m−1) increases, the averagemagnitude of each element of χ_(p) _(m) should decrease, if multipathray estimates are properly obtained. Given the size of R_(p) and{circumflex over (R)}_(p) _(m−1) , it is also known that the size ofχ_(p) _(m) . The indices y and z are used to denote the row number andthe column number, respectively, for each element in each of these threematrices.

To obtain initial estimates for the Doppler frequency offset {tilde over(f)}_(p) _(m) , the propagation delay {tilde over (τ)}_(p) _(m) , andthe carrier phase {tilde over (φ)}_(p) _(m) , the following operation isfirst executed:

$\begin{matrix}{\left( {y_{peak},z_{peak}} \right) = {\underset{y,z}{argmax}\left\{ \left\{ {{\chi \left( {y,z} \right)}} \right\}_{y = 0}^{N_{\tau}} \right\}_{z = 0}^{N_{f} - 1}}} & (49)\end{matrix}$

The following three equations are now used to obtain {tilde over(f)}_(p) _(m) , {tilde over (τ)}_(p) _(m) , and {tilde over (φ)}_(p)_(m) :

{tilde over (f)} _(p) _(m) =β_(y) _(peak)   (50)

{tilde over (τ)}_(p) _(m) =γ_(z) _(peak)   (51)

{tilde over (φ)}_(p) _(m) =∠χ(y _(peak) ,z _(peak))  (52)

where ∠ denotes the phase angle computation, which is the same as thefour quadrant arctangent operationtan⁻¹(imag(χ(y_(peak),z_(peak)))/real(χ(y_(peak),z_(peak)))).

Next a computation of an initial amplitude estimate is performed inblock 44. This process accepts {tilde over (f)}_(p) _(m) , {tilde over(τ)}_(p) _(m) , and {tilde over (φ)}_(p) _(m) as inputs and outputs aninitial estimate of the peak amplitude of the ray being considered,Ã_(p) _(m) . The initial amplitude is obtained through the use of knownlinear regression algorithms that have been found to sufficientlyestimate the peak amplitude of the ray of interest, m.

The search space is decomposed in block 46. This process accepts Ã_(p)_(m) , {tilde over (f)}_(p) _(m) , {tilde over (τ)}_(p) _(m) , and{tilde over (φ)}_(p) _(m) as inputs and outputs a vector containing awaveform design making use of the parameter estimates for ray m,{circumflex over (r)}_(p) _(m) , to the process used to compute thesearch space for the estimate as set forth in further detail below. Thefinal parameters associated with {circumflex over (r)}_(p) _(m) are thepeak amplitude estimate Ã_(p) _(m) , the frequency estimate {tilde over(f)}_(p) _(m) , the propagation delay estimate {tilde over (τ)}_(p) _(m), and the carrier phase estimate {tilde over (φ)}_(p) _(m) and areoutput to the ray database as also set forth in further detail below.

FIG. 3 illustrates a flowchart 48 outlining a process used to decomposethe search space. This flowchart 48 begins at block 50 and involves theexecution of several subprocesses. Generally, these subprocessesinclude: defining of the boundaries for the ranges of trial parametervalues, computation of the search grid for received data that will beused in the decomposition subprocess, construction of the code delayvector that will be used in designing the trial ray waveform,construction of the trial ray waveform, computation of the search gridresulting from addition of the trial ray waveform with previouslyobtained multipath ray estimates, comparison of the received search gridwith the search grid computed using estimated and trial parameters,computation of the error using the trial parameters as estimates of themultipath ray parameters of interest, and generation of new trialparameters.

Boundaries are defined for trial parameters in block 52. This subprocessaccepts Ã_(p) _(m) , {tilde over (f)}_(p) _(m) , {tilde over (τ)}_(p)_(m) , and {tilde over (φ)}_(p) _(m) as inputs and outputs boundaryvalues (both minimum and maximum) for the trial parameters of the raybeing considered. Parameter boundaries are output for the peak amplitude(A_(min) and A_(max)), the carrier frequency (f_(min) and f_(max)), thepropagation delay (τ_(min) and τ_(max)), and the carrier phase (φ_(min)and φ_(max)). Also made available from previous processes are r_(p) andN_(p).

This subprocess is used to simply define the ranges from which trialparameters for the peak amplitude A_(trial), the carrier frequencyf_(trial), the propagation delay τ_(trial), and the carrier phaseφ_(trial) in the decomposition process may be obtained. These ranges aredefined as follows:

A _(trial) ε[A _(min) ,A _(max)]  (53)

f _(trial) ε[f _(min) , f _(max)]  (54)

τ_(trial)ε[τ_(min),τ_(max)]  (55)

φ_(trial)ε[φ_(min),φ_(max)]  (56)

In this illustrated embodiment of the invention, these boundaries aredefined as follows: A_(min)=0 and A_(max)=max{|r_(p)|}. This range isset to simply dictate in generation of trial parameters that the peakamplitude of the ray waveform must not be outside the range of possiblevalues found in the received data vector. Because of the signal to noiseratio of GPS data recorded using a SiGe or TRIGR receiver, this range isappropriate, since peak noise amplitude values will almost certainlyexceed the peak amplitude of a GPS waveform. For some simulationcircumstances, however, A_(max) may need to be modified. The carrierfrequency trial parameter range boundaries are f_(min)={tilde over(f)}_(p) _(m) −500+Δf and f_(max)={tilde over (f)}_(p) _(m) +500−Δf. Asafety margin of Δf is used to bound estimates to within (but notincluding) a 1 kHz null-to-null bandwidth of a GPS L1 signal. The rangeof approximately −500 to 500 Hz around the initial estimate {tilde over(f)}_(p) _(m) is used because the absolute difference between two GPSwaveform frequency responses reach a peak at approximately 500 Hz. Thepropagation delay trial parameter range boundaries are τ_(min)={tildeover (τ)}_(p) _(m) −T_(c)+Δτ and τ_(max)={tilde over (τ)}_(p) _(m)+T_(c)−Δτ. A safety margin of Δτ is used to bound estimates to within(but not including) a chip duration of a GPS L1 signal. The carrierphase trial parameter range boundaries are φ_(min)=−2π and φ_(max)=2π.This range is chosen to ensure the true value of the carrier phase fallsin the range, and to account for the circular, rather than linear,nature of carrier phase values (a carrier phase of −π is equal to acarrier phase of π).

A received data search grid computation is then made in block 54. Thissubprocess accepts r_(p) and N_(p) as inputs and outputs a search gridto be used in the search space decomposition R_(s). Additional outputsof this subprocess are initial values for a trial iteration index g=0and a best error computation, E_(best)=∞. Graph 56 in FIG. 4Aillustrates the magnitude of the points and graph 58 in FIG. 4Billustrates the phase of the points in the search grid for simulatedreceived data R_(s).

Note that R_(s) is distinct from R_(p). This is a design choice that hasbeen made for this illustrated embodiment, as the initial estimate usedin the search space decomposition process is required to be relativelyaccurate, thus the use of a high resolution. The search grid R_(s) isused to obtain the estimate for the key parameters in an iterativeprocess. The smaller the resolution of the search grid, the greater thecomputational requirement. Therefore, construction of these two searchgrids is decoupled to ensure the simultaneous accuracy of initialestimates while still retaining an ability to perform the search spacedecomposition in a reasonable timeframe.

The search grid R_(s) used in a stochastic search process is computed asfollows:

R _(S)=grid(r _(p),τ_(window) ,f _(window),Δτ_(input) ,[Δf _(s)]_(N)_(fft) _(=N) _(p) )  (57)

where τ_(window)=2T_(c) (as before), f_(window)=10,000 kHz (as before),Δτ_(input)=T_(s) (as before), and the frequency resolution Δf_(s) issuch that the number of points in an FFT used to compute the search gridis equal to N_(p). Eq. 36 is used in the computation of Δf_(s). Thenumber of code delay values, thus the number of rows in the search grid,remains N_(τ). The number of frequency values, thus the number ofcolumns in the search grid, is denoted N_(f) _(s) .

Returning to FIG. 3, next a trial code delay vector is constructed inblock 60. This subprocess accepts R_(s), g, and E_(best) as inputs andoutputs a code delay vector x_(trial). This subprocess additionallyaccepts A_(trial), f_(trial), τ_(trial), and φ_(trial) in instanceswhere g>0. From this point forward, the following trial parameters aredefined as follows, for cases when g=0:

A _(trial) =Ã _(p) _(m)   (58)

f _(trial) ={tilde over (f)} _(p) _(m)   (59)

τ_(trial)={tilde over (τ)}_(p) _(m)   (60)

φ_(trial)={tilde over (φ)}_(p) _(m)   (61)

For cases when g>0, the initial estimates for the four parameters are nolonger used, making use only of the trial parameters generated in thesubprocess described in further detail below.

x_(trial) is computed using Eq. 40, where τ_(s) in Eq. 40 is as follows:τs=τ_(trial). N_(s) in Eq. 40 is set equal to N_(p).

Next a trial ray waveform vector is constructed in block 62. Thissubprocess accepts x_(trial) as an input and outputs a trial raywaveform vector r_(trial), as well as a cumulative discrete-timewaveform containing a summation of the trial ray waveform with allpreviously obtained multi path ray estimates, r_(grid). r_(trial) iscomputed using Eq. 23, where A_(p) _(m) =A_(trial), D_(p) _(m) =1, x_(p)_(m) =x_(trial), f_(p) _(m) =f_(trial), and φ_(p) _(m) =exp(jφ_(trial)).The trial carrier frequency vector f_(trial) is defined as follows:

f _(trial)=[exp(j2πf _(trial) ωT _(s))]_(ω=0) ^(N) ^(p) ⁻¹  (62)

r_(trial) must be considered for instances when the received andrecorded data is real-valued (r_(p)ε

) versus when the recorded data is complex-valued (r_(p)ε

). The use of Eq. 23 applies directly when r_(p)ε

. If r_(p)ε

, r_(trial) is expressed as follows:

r _(trial) =Re{r _(trial) },r _(p)ε

  (63)

The decision point “m=0” in decision block 64 is used to determinewhether or not the ray waveforms from previous estimates should be addedto the trial waveform r_(trial). If m=0 (“Yes” branch of decision block64), then it is assumed the LOS ray is being addressed. If another rayin the PIT is being addressed, then it is assumed a non-LOS ray is beingaddressed, and m>0. However, each ray estimate must be considered toreduce the difference between R_(p) and the estimate. Therefore, thewaveforms generated using previous ray estimates (resulting fromprevious iterations of the decomposition process illustrated in FIG. 2)must be added to the current ray waveform, in order to obtain anaccurate estimate of the current ray parameters. Therefore, if m>0 (“No”branch of decision block 64), the following operation is executed inblock 66, corresponding with the subprocess named for the followingequation, and addressed earlier:

$\begin{matrix}{{r_{grid} = {r_{trial} + {\sum\limits_{h = 0}^{m - 1}\; r_{p\; h}}}},{m > 0}} & (64)\end{matrix}$

If m=0, the following is true:

r _(grid) =r _(trial) ,m=0  (65)

A trial estimate search grid computation is performed after theconstruction of the trail ray waveform vector in block 62 if m=0 orafter the operation in block 66 if m≠0. This subprocess accepts r_(grid)as an input and outputs the search grid for the trial estimateR_(trial). Graph 70 in FIG. 5A illustrates the magnitude of the pointsand graph 72 in FIG. 5B illustrates the phase of the points in thesearch grid for the trial R_(trial).

The search grid is computed as follows:

R _(trial)=grid(r _(grid),τ_(window) ,f _(window),Δτ_(input) ,[Δf_(s)]_(N) _(fft) _(=N) _(p) )  (66)

Received and estimate search grids are then compared in block 74. Thissubprocess accepts R_(trial) as an input and outputs a comparison matrixχ_(trial). Graph 76 in FIG. 6A illustrates the magnitude of the pointsand graph 78 in FIG. 6B illustrates the phase of the points in thesearch grid for the difference between the received signal and thetrial, χ_(trial).

χ_(trial) is computed using the following equation:

χ_(trial) =R _(s) −R _(trial)  (67)

Returning again to FIG. 3, estimate Error is then computed in block 80.This subprocess accepts χ_(trial) as an input and outputs a ray estimateerror E_(trial) associated with the trial parameters. E_(trial) iscomputed using the following equation:

$\begin{matrix}{E_{trial} = {\frac{1}{N_{\tau}N_{f_{s}}}{\sum\limits_{y = 0}^{N_{\tau} - 1}\; {\sum\limits_{z = 0}^{N_{f_{s}} - 1}{{\chi_{trial}\left( {y,z} \right)}}^{2}}}}} & (68)\end{matrix}$

The decision point “E_(trial)<E_(max)” in decision block 82 is used todetermine whether or not the current error E_(trial) is superior toprevious error values. If this is the case (“Yes” branch of decisionblock 82), the parameters for the best estimate are asserted as thecurrent parameters in the subsequent subprocess in block 84 and denoted“E_(best)←E_(trial), r_(best)←r_(trial), A_(best)←A_(trial),τ_(best)←τ_(trial), f_(best)←f_(trial), φ_(best)←φ_(trial)”.

The iteration index g is incremented in the process that follows theseoperations in block 86.

The decision point “g<g_(max)” in decision block 88 is used to determinewhether further iterations of trial parameter generation and comparisonwith received data shall take place. If g<g_(max) (“Yes” branch ofdecision block 88), then further iterations occur at block 90. Otherwise(“No” branch of decision block 88), the best estimates are asserted asthe final parameters in block 92 and the process ends at block 94. Inthis illustrated embodiment, the number of iterations used is typicallyg_(max)=1000 when using simulated annealing to generate new trialparameter values, though g_(max)=1000 is an arbitrary value typicallyused for the number of iterations in generic simulated annealingliterature. If other methods are used for stochastic search andestimation in other embodiments, a convergence threshold may be usedinstead.

The process to assert the best estimates as the final estimates, asdiscussed earlier, is denoted “Â_(p) _(m) ←A_(best), {circumflex over(f)}_(p) _(m) ←f_(best), {circumflex over (τ)}_(p) _(m) ←τ_(best),{circumflex over (φ)}_(p) _(m) ←φ_(best)”. Note that final estimatevalues will always be output from the search space decompositionprocess, regardless of the quality of the estimate.

Graphs 96 a-96 b, 96 c, and 96 d in FIG. 7 illustrate a set of varioustrial parameters generated using the simulated annealing process, asdescribed below. Note that the use of reannealing between iterations 600to 700. Reannealing is a restarting of the annealing process to testwhether the global optimum is reached in earlier iterations. Graph 98 inFIG. 8 illustrates the error associated with each set of these trialparameters.

New trial parameter values are generated in block 90. This subprocessaccepts g as an input and outputs new values of A_(trial), f_(trial),τ_(trial), and φ_(trial). There are different means by which to computetrial parameters, depending on the stochastic search and optimizationmethod used. This illustrated embodiment utilizes simulated annealingfor selection of trial parameters. There may, however, be otherstochastic search and optimization techniques that may be employed formore accurate or efficient performance of the algorithm that may be usedin other embodiments of the invention.

Simulated annealing is made available to users practically in theMATLAB® global optimization toolbox, and has been employed in thisillustrated embodiment to serve as the computational backbone for trialparameter generation. Specific simulated annealing parameters associatedwith the generation of initial research findings will be outlined withthe results.

Returning now to FIG. 2, the estimate search space may now be computedin block 100. This process accepts {circumflex over (r)}_(p) _(m) as aninput and outputs an estimate search grid {circumflex over (R)}_(p) _(m). The input {circumflex over (r)}_(p) _(m) is output from the searchspace decomposition process outlined above and illustrated in FIG. 3.

The search grid is computed using the following equation:

$\begin{matrix}{{\hat{R}}_{p_{m}} = {{grid}\left( {{\sum\limits_{h = 0}^{m}\; {\hat{r}}_{p\; h}},\tau_{window},f_{window},{\Delta\tau}_{input},{\Delta \; f_{input}}} \right.}} & (69)\end{matrix}$

Note with this operation that the search grid obtained using the sum ofall ray waveforms must be computed.

After the computation of the estimate search space, a comparison ofreceived and estimate search spaces is performed in block 102. Thisprocess accepts {circumflex over (R)}_(p) _(m) as an input and outputs acomparison matrix χ_(p) _(m) . The equation used to compute χ_(p) _(m)is as follows:

χ_(p) _(m) =R _(p) _(m) −{circumflex over (R)} _(p) _(m)   (70)

This process is executed distinctly from the comparison subprocessdescribed in conjunction with block 74 of FIG. 3 as part of the searchspace decomposition process. While the result of the comparison processin block 74 of FIG. 3 would yield a similar result, the processdescribed above in conjunction with block 102 is different in that thesearch grids denoted χ_(p) _(m) , R_(p) _(m) , and {circumflex over(R)}_(p) _(m) may use a different propagation delay resolution Δτ and adifferent frequency resolution Δf from that used in block 74 of FIG. 3.

The estimate error is then computed in block 104. This process acceptsχ_(p) _(m) as an input and outputs the estimate error E_(p) _(m) . Graph106 in FIG. 9 illustrates the error associated with each ray estimate.The true number of multipath rays used in constructing graph 106 isfive, though the error reduction between the estimates for the fifth andsixth rays may lead one to believe that there are six multipath rays.This results simply as a function of an erroneous estimate thatcoincidentally significantly reduces the error.

The error E_(p) _(m) is computed as follows:

$\begin{matrix}{E_{p_{m}} = {\frac{1}{N_{\tau}N_{fs}}{\sum\limits_{y = 0}^{N_{\tau} - 1}\; {\sum\limits_{z = 0}^{N_{fs} - 1}\; {{\chi_{p_{m}}\left( {y,z} \right)}}^{2}}}}} & (71)\end{matrix}$

The decision point “m=0” in decision block 108 is used to determinewhether the first, or LOS, ray is being considered. If so (“Yes” branchof decision block 108), the estimate parameters output from the searchspace decomposition are taken as is to be the estimate parameters forthe LOS ray. If m>0 (“No” branch of decision block 108), advance is madeto another decision point, denoted “E_(p) _(m) −E_(p) _(m−1) <ε” indecision block 110. This decision is made to determine if the currentray estimate improves the global estimate that makes use of all rayparameters obtained using the decomposition algorithm. Therefore, if theestimate associated with ray m improves the estimate error over theestimate associated with ray m−1 by at least the amount specified by theacceptance threshold ε (“Yes” branch of decision block 110), then theestimate is declared to be a valid parameterization of the current ray.If the estimate does not improve the error by at least ε (“No” branch ofdecision block 110), then the estimate parameters are not considered tobe valid, the estimate parameters are discarded, and other regionalmaxima are considered in block 112. The acceptance threshold used in theexample illustrated in graph 106 in FIG. 9 is ε=0.

If the estimate associated with ray m improves the estimate error overthe estimate associated with ray m−1 by at least the amount specified bythe acceptance threshold ε (“Yes” branch of decision block 110), theparameters are output to the ray database in block 114. This processaccepts Ã_(p) _(m) , {tilde over (f)}_(p) _(m) , {tilde over (τ)}_(p)_(m) , and {tilde over (φ)}_(p) _(m) as inputs and does not outputanything. The ray database may be used to simply store estimateparameters for use in inter-PIT algorithms. The decomposition algorithmis utilized to extract these parameters, eventually yielding thisdatabase. Using this database, multipath rays present in more than onePIT may be considered.

After this process is completed, the ray index m is incremented in block116, and then another iteration of the decomposition algorithm isperformed beginning again at block 42.

If the estimate above does not improve the error by at least ε (“No”branch of decision block 110), alternate regional maxima may bedetermined in block 112. This process accepts χ_(p) _(m) ⁻¹ as an input.The absolute values of each of the elements in χ_(p) _(m) are evaluatedto determine the locations of regional maxima. The MATLAB®“imregionalmax” command may be used to determine these locations insoftware, and the 8-neighborhood regional maxima are determined. Ifthere is a regional maximum within χ_(p) _(m) that has not yet beenconsidered for estimation (“Yes” branch of decision block 118), then aninitial estimation of these alternate peak parameters takes place inblock 120. Otherwise (“No” branch of decision block 118), the PIT indexis incremented in block 122 and the next PIT is decomposed. If p>p_(max)(“Yes” branch of decision block 124), there are no more PIT-sizedvectors of received data to consider, and the algorithm is completed atblock 126. Otherwise (“No” branch of decision block 124), the processcontinues with the next PIT-sized vector at block 38.

If there is a regional maximum within χ_(p) _(m) that has not yet beenconsidered for estimation (“Yes” branch of decision block 118), then aninitial estimation of these alternate peak parameters takes place inblock 120. This process accepts χ_(p) _(m−1) as an input and outputs{tilde over (f)}_(p) _(m) , {tilde over (τ)}_(p) _(m) , and {tilde over(φ)}_(p) _(m) . The choices for the parameters {tilde over (f)}_(p) _(m), {tilde over (τ)}_(p) _(m) , and {tilde over (φ)}_(p) _(m) correspondwith the frequency, the propagation delay, and the carrier phase for theelement in χ_(p) _(m−1) with the highest magnitude that is a regionalmaximum and has not yet been considered. These parameters are then usedto estimate the initial amplitude and the process continues at block 44.

Simulations are conducted to determine the accuracy of the abovedescribed signal decomposition and parameterization algorithm inembodiments of the invention for obtaining estimates for ray parameters,in an attempt to quantify the performance of the algorithm. In thesimulations, parameter estimates are compared with true parametervalues. The following metrics are used to quantify performance:

-   -   The bias and standard deviation in amplitude estimates relative        to the true amplitude values    -   The bias and standard deviation in Doppler frequency offset        estimates relative to the true Doppler frequency offset values    -   The bias and standard deviation in propagation path length        offset estimates relative to the true propagation path length        values    -   The bias and standard deviation in carrier phase estimates        relative to the true carrier phase values

Simulations are conducted using the GPS waveform model. To determine theeffectiveness of the signal decomposition and parameterization algorithmin an environment with a varying C/N₀, the simulated GPS signal and thesimulated noise must be scaled before summing the signal and the noisetogether. This scaling is performed in accordance with theory found inliterature, for example Misra et al., “Global Positioning System:Signals, Measurements, and Performance,” 2d Ed., Lincoln, Mass.,Ganga-Jamuna Press, 2006. Only PRN 26 is considered in simulations,omitting signals from other PRNs normally available in received data.The choice of PRN 26 is entirely arbitrary. The presence of an automaticgain controller is simulated in order to consider data that takesadvantage of the full dynamic range afforded by the number of bits usedin generating a data set without clipping signal peaks.

Simulations are conducted for various values of C/N₀: 40, 45, and 50dB-Hz, as well as simulations where no noise is present. The number ofquantization levels is varied in the simulations as well: 1-bit, 4-bit,8-bit, and 16-bit fixed, as well as double-precision floating point. Foreach case of a quantization level coupled with a value of C/N₀, 1000integration periods using a PIT of 1 msec are generated where onediscrete-time ray vector, presumed in processing to be line-of-sight,from PRN 26 is added to white Gaussian noise. The simulated intermediatefrequency is −38.4 kHz and the simulated sampling frequency is 8.1838MHz, in line with parameters that would be used for a SiGe GPS frontend. The path length offset of the ray within each integration period isfixed, varies between integration periods, and is selected arbitrarilyusing a uniform random distribution ranging from −cT_(c) to cT_(c). TheDoppler frequency offset of the ray in each integration period is fixed,varies between integration periods, and is selected arbitrarily using auniform random distribution ranging from −8 kHz to 8 kHz. The initialcarrier phase of the ray in each integration period is fixed, variesbetween integration periods, and is selected arbitrarily using a uniformrandom distribution ranging from −180° to 180°. For every ray estimate,1000 simulated annealing iterations are used to obtain the estimate,i.e., g_(max)=1000.

FIGS. 10-13 illustrate the normalized amplitude estimate errordistribution (graph 128 in FIG. 10), the carrier frequency offsetestimate error distribution (graph 130 in FIG. 11), the propagation pathlength offset error distribution (graph 132 in FIG. 12), and the carrierphase error distribution (graph 134 in FIG. 13), respectively, for thesimulation conducted using double precision floating point numbers andC/N₀=50 dB-Hz. These figures are provided to serve as examples oftypical error distribution shapes found in the histograms resulting fromthe simulations set out above.

The table in FIG. 14 illustrates a comparison of the mean amplitudeerror (normalized) of simulated rays for various quantization levels andvarious values of C/N₀. The table in FIG. 15 illustrates a comparison ofthe standard deviation of the amplitude error (normalized).Normalization in this context is defined by the following formula:(estimated value−true value)/true value. The table in FIG. 16illustrates a comparison of the mean frequency offset error (Hz) ofsimulated rays for various quantization levels and various values ofC/N₀. The table in FIG. 17 illustrates a comparison of the standarddeviation of the frequency offset error (Hz). The table in FIG. 18illustrates a comparison of the mean path length offset error (m) ofsimulated rays for various quantization levels and various values ofC/N₀. The table in FIG. 19 illustrates a comparison of the standarddeviation of the path length offset error (m). The table in FIG. 20illustrates a comparison of the mean carrier phase error (deg) ofsimulated rays for various quantization levels and various values ofC/N₀. The table in FIG. 21 illustrates a comparison of the standarddeviation of the carrier phase error (deg).

The results provided assist in establishing an understanding of theperformance of the algorithm set out above. The results indicate thatthe algorithm is indeed able to obtain accurate estimates for each ofthe four key parameters used to define a ray in a setting with a higherC/N₀. Higher dynamic ranges demonstrate superior potential forestimation accuracy. Both of these results are to be expected, as bothof these results favor the scenario where less noise is present.

The amplitude error distribution (Graph 128), the carrier frequencyoffset error distribution (Graph 130), and the carrier phase errordistribution (Graph 134) all appear to share similarly shapeddistribution function (with the exception of the outlier cases where theamplitude is improperly, as observed in Graph 128). Neglecting thisrelatively small percentage of amplitude errors, all three of thesedistributions appear to be triangular in nature, though a Gaussiandistribution may perhaps adequately describe these error distributions.The propagation path length offset error distribution (Graph 132)appears, on the other hand, to have a uniform distribution with a biasof approximately −14 m among simulations yielding reasonable results.Furthermore, the distribution appears in histograms to range from −35 to+5 m. The path length that would be observed between rays that arereceived exactly one sampling period T_(s) apart from each other, usingthe sampling frequency used in these simulations, is 36.63 m. This isrelatively close to the approximately 40 m size of the uniformdistribution of results. It has been observed that path length offsetestimates using the algorithm's simulated annealing function to estimateparameters are very similar to the initial estimates obtained toinitialize the simulated annealing function. Therefore, use of thesignal decomposition and parameterization algorithm does not appear toprovide better path length offset estimates than what would be obtainedusing the coarse estimates presented by the search grid for initialestimation and obtained using Eq. 51.

Flowchart 136 in FIG. 22 illustrates a top-level flow of an alternateembodiment of the intra-PIT signal decomposition algorithm. Flowchart136 corresponds directly with the process flowchart 138 illustrated inFIGS. 23A and 23B. The process begins at block 140. Similar to theprevious illustrated embodiment, RF front end functions (such asfrequency downconversion and signal sampling) are performed at block142. The signal is then aligned such that data is processed by adecomposition algorithm in integer multiples of the duration of the PRNcode sequence present in GPS L1 C/A-coded signals, for example, ininteger multiples of 1 msec. As with the previous illustratedembodiment, the decomposition algorithm may only to process wholeperiods of the PRN code sequence. After alignment, received data that issized to be of PIT-sized duration is input for search grid computationand initially ray parameter estimation in block 146. Amplitude, Dopplerfrequency, relative propagation delay, and carrier phase are thencoarsely estimated. These coarse parameter estimates may then be used toinitiate fine estimation of these four parameters for a ray waveformusing a stochastic search and optimization technique (simulatedannealing) in some embodiments (also in block 146), and a portion of thesearch grid associated with the estimate parameters may then be removedfrom the search grid.

The remainder of the search grid may then be evaluated to obtain astopping criteria statistic in block 148. The stopping criteriastatistic may be obtained in some embodiments by comparing a peak powerof the search grid with a noise power present in the search grid, thoughother criteria in other embodiments may also be used. The stoppingcriteria statistic is then evaluated to determine if the stoppingcriteria is satisfied in block 150. If the criteria is satisfied (“Yes”branch of decision block 150), then all rays present in the search gridare treated as having been estimated by the algorithm (block 152), andthe algorithm advances to the next PIT for decomposition in block 154.If the criteria is not satisfied (“No” branch of decision block 150),then the search grid is evaluated to compute an estimate error in block156, and the algorithm processor decides if the estimate error can befurther reduced in block 158. If so (“Yes” branch of decision block158), decomposition of the current PIT continues at block 160, and amatrix difference between the received and estimate search grids iscomputed in block 162 for use in subsequent coarse ray parameterestimation. This process continues until there are no more rays toestimate. If error cannot be reduced (“No” branch of decision block158), then all rays present in the search grid are treated as havingbeen estimated by the algorithm in block 152, and the algorithm advancesto the next PIT for decomposition in block 154. This loop is repeateduntil all data is processed and parameterized.

The previous illustrated embodiment utilizes hard estimate decisions onray waveform parameter estimates. The current illustrated embodimentutilizes soft ray waveform parameter estimate decisions that areobtained simultaneously through the parallel estimation of multiple raywaveform parameters, as part of the stochastic search block 146 in FIG.22. This approach to ray waveform parameter estimation may be referredto as “multiray” estimation.

The description of the algorithm illustrated in flowchart 138 in FIGS.23A and 23B utilizes the continuous time-domain signal model of Eq. 1set out above. The process begins at block 164. Received data isdownconverted at block 166. The downconversion may be executed in RFfront end hardware of a software receiver similar to the previousembodiment illustrated above. After the downconversion process, data issampled at block 168. PIT index value, p, is initialized to zero inblock 170. After sampling, the Code-aligned data may be buffered atblock 172. ε is initialized to zero at block 174. A search grid iscomputed in block 176. An initial search of estimated grid peaklocations is performed in block 178 and initial parameter estimation isperformed in block 180. This series of steps utilize the samemethodology as set out in the previous illustrated embodiment.

FIG. 24 illustrates a flowchart 184 outlining an alternate process usedto decompose the search space of block 182 in flowchart 136. Thisalternate process accepts either Ã_(p) _(ε) , {tilde over (f)}_(p) _(ε), {tilde over (τ)}_(p) _(ε) , and {tilde over (φ)}_(p) _(ε) (in caseswhen ε=0) or Ã_(p) _(ε) , {tilde over (f)}_(p) _(ε) , {tilde over(τ)}_(p) _(ε) , and {tilde over (φ)}_(p) _(ε) (in cases when ε>0).Output from this process are the final ray ensemble parameter estimatesassociated with {circumflex over (r)}_(p) _(ε) : the peak amplitudeestimates Â_(p) _(ε) , the frequency estimate {circumflex over (f)}_(p)_(ε) , the propagation delay estimates {circumflex over (τ)}_(p) _(ε) ,and the carrier phase estimates {circumflex over (φ)}_(p) _(ε) .Estimates are obtained for each of the ε+1 rays in ensemble ε.

Ã_(p) _(ε) , {tilde over (f)}_(p) _(ε) , {tilde over (τ)}_(p) _(ε) ,{tilde over (φ)}_(p) _(ε) , Â_(p) _(ε) , {circumflex over (f)}_(p) _(ε), {circumflex over (τ)}_(p) _(ε) , and {circumflex over (φ)}_(p) _(ε)are defined as follows:

Ã _(p) _(ε) =[Â _(p) _(ε−1) ,Ã _(p) _(ε) ]  (72)

{tilde over (f)} _(p) _(ε) =[{circumflex over (f)} _(p) _(ε−1) ,{tildeover (f)} _(p) _(ε) ]  (73)

{tilde over (τ)}_(p) _(ε) =[{circumflex over (τ)}_(p) _(ε−1) ,{tildeover (τ)}_(p) _(ε) ]  (74)

{tilde over (φ)}_(p) _(ε) =[{circumflex over (φ)}_(p) _(ε−1) ,{tildeover (φ)}_(p) _(ε) ]  (75)

where

Â _(p) _(ε) =[Â _(p) _(ε) ]_(m=0) ^(ε)  (76)

{circumflex over (f)} _(p) _(ε) =[{circumflex over (f)} _(p) _(ε)]_(m=0) ^(ε)  (77)

{circumflex over (τ)}_(p) _(ε) =[{circumflex over (τ)}_(p) _(ε) ]_(m=0)^(ε)  (78)

{circumflex over (φ)}_(p) _(ε) =[{circumflex over (φ)}_(p) _(ε) ]_(m=0)^(ε)  (79)

The decomposition process begins at block 186. Boundaries are definedfor trial parameters in block 188. This subprocess accepts either Ã_(p)_(ε) , {tilde over (f)}_(p) _(ε) , {tilde over (τ)}_(p) _(ε) , and{tilde over (φ)}_(p) _(ε) (in cases when ε=0) or Ã_(p) _(ε) , {tildeover (f)}_(p) _(ε) , {tilde over (τ)}_(p) _(ε) , and {tilde over(φ)}_(p) _(ε) (in cases when ε>0). Output from this sub-process areboundary values (both minimum and maximum) for the trial parameters usedto generate trial ray waveforms being considered. Parameter boundariesare output for the peak amplitude (A_(min) and A_(max)), carrierfrequency (f_(min) and f_(max)), propagation delay (τ_(min) andτ_(max)), and the carrier phase (φ_(min) and φ_(max)). Also madeavailable is N_(p). This subprocess is used to simply define the rangesfrom which trial parameters for the peak amplitude A_(trial), thecarrier frequency f_(trial), the propagation delay τ_(trial), and thecarrier phase φ_(trial) in the decomposition process may be obtained.These ranges are defined as follows:

A _(trial) ε[A _(min) ,A _(max)]  (80)

f _(trial) ε[f _(min) ,f _(max)]  (81)

τ_(trial)ε[τ_(min),τ_(max)]  (82)

φ_(trial)ε[φ_(min),φ_(max)]  (83)

where trial parameter vectors A_(trial), f_(trial), τ_(trial), andφ_(trial) generated using simulated annealing are expressed as follows:

A _(trial) =[A _(trial) _(m) ]_(m=0) ^(ε)  (84)

f _(trial) =[f _(trial) _(m) ]_(m=0) ^(ε)  (85)

τ_(trial)=[τ_(trial) _(m) ]_(m=0) ^(ε)  (86)

φ_(trial)=[φ_(trial) _(m) ]_(m=0) ^(ε)  (87)

In this embodiment, boundaries are defined as follows: A_(min)=0 andA_(max)=max{|r_(p)|}. This particular range was selected to dictate ingeneration of trial parameters that a peak amplitude of a ray waveformmay not be outside a range of possible values found in the received datavector. Due to a signal to noise ratio of GPS data recorded using asoftware receiver, this range is appropriate since peak noise amplitudevalues will almost certainly exceed the peak amplitude of a GPSwaveform. For simulation circumstances in some embodiments, however,A_(max) may need to be modified. Carrier frequency parameter rangeboundaries and propagation delay parameter range boundaries are definedsimilar to the previously illustrated embodiment. These boundaries mayalso be employed for this subprocess. The carrier frequency trailparameter range boundaries are f_(min)={tilde over (f)}_(p) ₀−1.1/T_(PIT) and f_(max)={tilde over (f)}_(p) ₀ +1.1/T_(PIT) in caseswhen ε=0 or f_(min)={circumflex over (f)}_(p) ₀ −1.1/T_(PIT) andf_(max)={circumflex over (f)}_(p) ₀ +1.1/T_(PIT) in cases when ε>0. Incases when ε>0, f_(min) and f_(max) are scalars

(f̂_(p₀) = f̂_(p_(0₀)))

since only one ray waveform is parameterized in the estimate waveformgenerated from the initial ensemble in block 180 of flowchart. Thus,this highest ray index for an ensemble is equal to the ensemble indexitself. The propagation delay trial parameter range boundaries areτ_(min)={tilde over (τ)}_(p) ₀ −1.1/T_(C) and τ_(max)={tilde over(τ)}_(p) ₀ +1.1/T_(C) in cases when ε=0 or τ_(min)={circumflex over(τ)}_(p) ₀ −1.1/T_(C) and τ_(max)={circumflex over (τ)}_(p) ₀ +1.1/T_(C)in cases when ε>0. In cases when ε>0, τ_(min) and τ_(max) are scalarssince

τ̂_(p₀) = τ̂_(p_(0₀)).

The carrier phase trial parameter range boundaries are φ_(min)=−2π andφ_(max)=2π. This range is chosen to ensure the true value of the carrierphase falls in the range, and to account for the circular, rather thanlinear, nature of carrier phase values, though other embodiments may useother ranges.

A received data search grid computation is then made in block 190. Thissubprocess accepts r_(p) and N_(p) as inputs and outputs a search gridto be used in the search space decomposition R_(s). Note that R_(s) isdistinct from R_(p). This is a design choice that has been made for thisembodiment, as an initial estimate used in the search spacedecomposition process is required to be relatively accurate, thus theuse of tight grid spacings, though other choices may be made for otherembodiments. The search grid R_(s) is used to obtain the estimate forthe key parameters in an iterative process. The tighter the grid spacingof the search grid, the greater the computational requirement.Therefore, construction of these two search grids is decoupled to ensurethe simultaneous accuracy of initial estimates while still retaining theability to perform the search space decomposition in a reasonabletimeframe.

Search grid R_(s) used in the stochastic search process is computed asfollows:

$\begin{matrix}{R_{s} = {{grid}\left( {{r_{p} \cdot \tau_{window}},f_{window},{\Delta\tau}_{input},\left\lbrack {\Delta \; f_{input}} \right\rbrack_{N_{fft} = N_{f_{input}}}} \right)}} & (87)\end{matrix}$

where τ_(window)=2T_(C), f_(window)=10 kHz, Δτ_(input)=T_(s), and thefrequency grid spacing Δf_(s) is such that the number of points in theFFT used to compute the search grid is equal to N_(f) _(input) . Theequation for N_(f) _(input) is as follows:

$\begin{matrix}{{N_{f_{input}} = 2^{\lambda}},{\lambda = \left\lceil {\log_{2}\left( {\max \left( {\frac{f_{s}}{N_{p}},\frac{2}{3N_{p}T_{s}}} \right)} \right)} \right\rceil}} & (88)\end{matrix}$

Note that Eq. 36 above is used in the computation of Δf_(s). The numberof code delay values, thus the number of rows in the search grid, isdenoted N_(τ) _(s) . The number of frequency values, thus the number ofcolumns in the search grid, is denoted N_(f) _(s) . After the searchgrid is computed, trial iteration index g is initialized to zero andbest error computation E_(best) is initialized to ∞ in block 192.

A trial estimate search grid is computed in block 194. This subprocessaccepts A_(trial), f_(trial), τ_(trial), and φ_(trial) as inputs andoutputs a search grid for the trial estimate R_(trial).

In the previously illustrated embodiment, the search grid is generatedthrough the computation of discrete time waveform r_(trial) and thedirect transformation of r_(trial) to the search grid R_(trial). Thismethod is more computationally intensive than the present embodiment,which utilizes an approximation to R_(trial) that follows:

$\begin{matrix}{R_{trial} \approx {{grid}\left( {r_{trial},\tau_{window},f_{window},{\Delta\tau}_{input},{\Delta \; f_{input}}} \right)} \approx} & \; \\\begin{bmatrix}{R_{trial}\left( {\beta_{0},\gamma_{0}} \right)} & {R_{trial}\left( {\beta_{0},\gamma_{1}} \right)} & \ldots & {R_{trial}\left( {\beta_{0},\gamma_{N_{f} - 1}} \right)} \\{R_{trial}\left( {\beta_{1},\gamma_{0}} \right)} & {R_{trial}\left( {\beta_{1},\gamma_{1}} \right)} & \; & {R_{trial}\left( {\beta_{0},\gamma_{N_{f} - 1}} \right)} \\\vdots & \; & \ddots & \vdots \\{R_{trial}\left( {\beta_{N_{\tau} - 1},\gamma_{0}} \right)} & {R_{trial}\left( {\beta_{N_{\tau} - 1},\gamma_{1}} \right)} & \cdots & {R_{trial}\left( {\beta_{N_{\tau} - 1},\gamma_{N_{f} - 1}} \right)}\end{bmatrix} & (89) \\{\mspace{79mu} {where}} & \; \\{{R_{trial}\left( {\beta_{y},\gamma_{z}} \right)} = {\sum\limits_{m = 0}^{ɛ}\; \left\lbrack {A_{{trial}_{m}}\sin \; {c\left( {m,z} \right)}{\varrho \left( {m,y} \right)}{\exp \left( {j\left( {{{\pi \left( {f_{{trial}_{m}} - \gamma_{z}} \right)}T_{PIT}} + \varphi_{{trial}_{m}}} \right)} \right)}} \right\rbrack}} & (90) \\{\mspace{79mu} {{\sin \; {c\left( {m,z} \right)}} = \begin{Bmatrix}{\frac{\sin \left( {{\pi \left( {f_{{trial}_{m}} - \gamma_{z}} \right)}T_{PIT}} \right)}{{\pi \left( {f_{{trial}_{m}} - \gamma_{z}} \right)}T_{PIT}},{{f_{{trial}_{m}} - \gamma_{z}} \neq 0}} \\{1,{{f_{{trial}_{m}} - \gamma_{z}} = 0}}\end{Bmatrix}}} & (91) \\{\mspace{79mu} {{\varrho \left( {m,y} \right)} \approx \begin{Bmatrix}{{1 - {\frac{\tau_{{trial}_{m}} - \beta_{y}}{T_{c}}}},{{{\tau_{{trial}_{m}} - \beta_{y}}} \leq T_{c}}} \\{0,{{{\tau_{{trial}_{m}} - \beta_{y}}} > T_{c}}}\end{Bmatrix}}} & (92)\end{matrix}$

The same grid parameters used to define the received search grid R_(s)in Eq. 87 are used to define the trial estimate search grid R_(trial) inEq. 89.

Received and estimate search grids are compared in block 196. Thissubprocess accepts R_(s) and R_(trial) as inputs and outputs acomparison matrix χ_(trial). χ_(trial) is computed by subtractingR_(trial) from R_(s). After the comparison, an error estimate iscomputed in block 198. This subprocess accepts χ_(trial) as an input andoutputs a trial estimate error E_(trial) associated with the trialparameters. E_(trial) is computed using the following equation:

$\begin{matrix}{E_{trial} = \sqrt{\frac{1}{N_{\tau_{s}}N_{f_{s}}}{\sum\limits_{y = 0}^{N_{\tau_{s}} - 1}\; {\sum\limits_{z = 0}^{N_{f_{s}} - 1}{{\chi_{trial}\left( {y,z} \right)}}^{2}}}}} & (93)\end{matrix}$

Decision point at block 200 is used to determine whether or not thecurrent error E_(trial) is superior to previous error values. If this isthe case (“Yes” branch of decision block 200), the parameters for thebest estimate are asserted as the current parameters in the subsequentsubprocesses at block 202. Regardless of whether or notE_(trial)<E_(best), the iteration index g is then incremented in block204. Decision point at block 206 is then used to determine whetherfurther iterations of trial parameter generation and comparison withreceived data shall take place. If g is less than g_(max) (“Yes” branchof decision block 206), then further iterations occur, starting with thegeneration of new trial parameter values at block 208. Otherwise (“No”branch of decision block 208), the best estimates are asserted as thefinal parameters in block 210. In this illustrated embodiment, thenumber of iterations (g_(max)) is set to 1000 when using simulatedannealing to generate new trial parameter values, though other numbersfor the maximum iterations may also be used. Upon assertion of finalparameter values in block 210, the search space decomposition processcompletes at block 212.

Returning now to flowchart 136 in FIGS. 23A and 23B, the estimate searchspace is computed in block 214. This subprocess accepts Â_(p) _(ε) ,{circumflex over (f)}_(p) _(ε) , {circumflex over (τ)}_(p) _(ε) , and{circumflex over (φ)}_(p) _(ε) as inputs and outputs an estimate searchgrid {circumflex over (R)}_(p) _(ε) . To compute {circumflex over(R)}_(p) _(ε) , an estimate discrete time-domain waveform {circumflexover (r)}_(p) _(ε) must first be computed using the following equation:

$\begin{matrix}{\mspace{79mu} {{{\hat{r}}_{p_{ɛ}} = {\sum\limits_{m = 0}^{ɛ}\; {\hat{r}}_{p_{ɛ_{m}}}}}\mspace{79mu} {where}}} & (94) \\{{\hat{r}}_{p_{ɛ_{m}}} = \left\{ \begin{matrix}{{{\hat{A}}_{p_{ɛ_{m}}}{{\hat{\chi}}_{p_{ɛ_{m}}} \cdot \left\lbrack {\exp \left( {j\left( {{2{\pi \left( {f_{IF} + {\hat{f}}_{p_{ɛ_{m}}}} \right)}T_{s}w} + {\hat{\varphi}}_{p_{ɛ_{m}}}} \right)} \right)} \right\rbrack_{w = 0}^{N_{p} - 1}}},} & {r_{p} \in {\mathbb{C}}} \\{{{Re}\left\{ {{\hat{A}}_{p_{ɛ_{m}}}{{\hat{\chi}}_{p_{ɛ_{m}}} \cdot \left\lbrack {\exp \left( {j\left( {{2{\pi \left( {f_{IF} + {\hat{f}}_{p_{ɛ_{m}}}} \right)}T_{s}w} + {\hat{\varphi}}_{p_{ɛ_{m}}}} \right)} \right)} \right\rbrack_{w = 0}^{N_{p - 1}}}} \right\}},} & {r_{p} \in {\mathbb{R}}}\end{matrix} \right.} & (95)\end{matrix}$

and the computation of

χ̂_(p_(ɛ_(m)))

makes use of Eq. 40 as follows:

$\begin{matrix}{{\hat{\chi}}_{p_{ɛ_{m}}} = {\chi_{s}\left( {\hat{\tau}}_{p_{ɛ_{m}}} \right)}} & (96)\end{matrix}$

To compute the search grid {circumflex over (R)}_(p) _(ε) , thefollowing equation is used:

{circumflex over (R)} _(p) _(ε) =grid({circumflex over (r)} _(p) _(ε),τ_(window) ,f _(window),Δτ_(input) ,Δf _(input))  (97)

Received and estimated search spaces are then compared in block 216.This subprocess accepts R_(p) and {circumflex over (R)}_(p) _(ε) asinputs and outputs the comparison matrix {circumflex over (χ)}_(p) _(ε)which is calculated by subtracting {circumflex over (R)}_(p) _(ε) from{circumflex over (R)}_(p) _(ε) . An error estimate is then computed inblock 218. This subprocess accepts χ_(p) _(ε) as an input and outputs anestimate error E_(p) _(ε) , which is computed as follows:

$\begin{matrix}{E_{p_{ɛ}} = \sqrt{\frac{1}{N_{\tau}N_{f}}{\sum\limits_{y = 0}^{N_{\tau} - 1}\; {\sum\limits_{z = 0}^{N_{f} - 1}\; {{\chi_{p_{ɛ}}\left( {y,z} \right)}}^{2}}}}} & (98)\end{matrix}$

Decision point in block 220 is used to determine whether the first, orLOS, ray is being considered. If ε=0 (“Yes” branch of decision block220), then the estimate parameters output from the search spacedecomposition are taken as is to be the estimate parameters for the LOSray in block 222. If ε>0 (“No” branch of decision block 220), then thealgorithm advances to another decision point at block 224. This decisionis made to determine if the current ensemble of ray parameter estimatesimproves the global estimate that makes use of all ray parameters withinε obtained using the decomposition algorithm. Therefore, if theparameter estimates associated with ε improve or remain the same inerror relative to the estimates associated with ensemble ε−1 (“Yes”branch of decision block 224), then the ensemble is declared to be avalid parameterization, and the ray parameters are output to the raydatabase in block 222. If the ensemble does not reduce or match theerror from the previous ensemble (“No” branch of decision block 224),then the ensemble is not considered to be valid, the ensemble isdiscarded, and other regional maxima are considered in blocks 226, 228,and 230 similar to the previous illustrated embodiment.

After parameters have been output to the ray database in block 222, astopping criteria statistic is computed in block 232. This subprocessaccepts χ_(p) _(ε) as an input and outputs stopping criteria statisticF_(p) _(ε) . F_(p) _(ε) may be described as a peak signal contentmagnitude to noise power comparison parameter. F_(p) _(ε) is computedusing a cumulative distribution function (cdf) of a Rayleighdistribution as follows:

$\begin{matrix}{{F_{p_{ɛ}} = {1 - {\exp \left( {- \frac{\xi_{p_{ɛ}}^{2}}{2\sigma_{x}^{2}}} \right)}}},{\xi \geq 0.}} & (99)\end{matrix}$

To obtain ξ_(p) _(ε) , Eq. 49 is used to determine the initial Dopplerfrequency propagation delay parameters, where y_(min) and y_(max) bindthe indices considered in χ_(p) _(ε) to those corresponding with thefollowing condition:

β_(y_(min)) ≤ τ̂_(p_(0₀)) ≤ β_(y_(max)).

z_(min) and z_(max) bind the indices considered in establishing ξ_(p)_(ε) to those corresponding with the following condition:

γ_(z_(min)) ≤ f̂_(p_(0₀)) ≤ γ_(z_(max)).

ξ_(p) _(ε) is then obtained using the following equation:

ξ_(p) _(ε) =|χ_(p) _(ε) (y _(peak) ,z _(peak))|  (100)

The carrier frequency trial parameter range boundaries aref_(min)={circumflex over (f)}_(p) ₀ −1.1/T_(PIT) and f_(max)={circumflexover (f)}_(p) ₀ +1.1/T_(PIT). The propagation delay trial parameterrange boundaries are τ_(min)={circumflex over (τ)}_(p) ₀ −(T_(c)+9×10⁻⁷)and τ_(max)={circumflex over (τ)}_(p) ₀ +(T_(c)+9×10⁻⁷). Multipathsignal content is expected to be delayed at most by 1.8×10⁻⁶ sec of theLOS ray, thus the use of this value in setting the boundary conditionsfor the propagation delay.

To obtain σ_(x), it is assumed that there is no signal content presentoutside the boundary conditions outlined in the previous paragraph. Thecomplex values of all the elements in χ_(p) _(ε) from grid points withindices corresponding with Doppler frequency or propagation delay valuesoutside these boundaries are concatenated into a single vector. Thestandard deviation σ_(x) of either the real or imaginary components isthen computed. Given the assumption that this noise content originatesfrom circular complex AWGN, the standard deviation of the realcomponents should be equal to the standard deviation of the imaginarycomponents.

Upon computation of F_(p) _(ε) , a stopping criteria is then appliedusing the following decision criteria in block 234:

$\begin{matrix}{F_{p_{ɛ}}\underset{Halt}{\overset{Continue}{\lessgtr}}\mathrm{\Upsilon}} & (101)\end{matrix}$

where

is the stopping criteria threshold. The stopping criteria threshold

may be determined through simulation or other methods in otherembodiments. If F_(p) _(ε) ≧

(“No” branch of decision block 234), then the ensemble index ε isincremented in block 236 and the next ensemble for the same PIT isconsidered for decomposition starting at block 178. If F_(p) _(ε) <

(“Yes” branch of decision block 234), decomposition of the current PITis halted and the PIT index is incremented in block 238. If there are nomore integration periods remaining to be processed (“Yes” branch ofdecision block 240), then the algorithm concludes at block 242.Otherwise (“No” branch of decision block 240), the next integrationperiod is process starting at block 172.

Embodiments of the signal decomposition and parameterization algorithmmay be implemented in any number of hardware/software configurations.FIG. 25 is a diagrammatic illustration of an exemplary hardware andsoftware environment for an apparatus 246 configured to execute thesignal decomposition and parameterization algorithm consistent withembodiments of the invention. Apparatus 246, in specific embodiments,may be a general purpose processor, application specific integratedcircuit (ASIC), field programmable gate array (FPGA), computer, computersystem, computing device, server, disk array, or programmable devicesuch as a multi-user computer, a single-user computer, a handheldcomputing device, a networked device (including a computer in a clusterconfiguration), a mobile telecommunications device, a video game console(or other gaming system), etc. Apparatus 246 may be referred to as“computing apparatus,” but will be referred to herein as a “generalpurpose processor.”

The general purpose processor 246 includes at least one centralprocessing unit (“CPU”) 248 coupled to a memory 250. Each CPU 248 istypically implemented in hardware using circuit logic disposed on one ormore physical integrated circuit devices or chips. Each CPU 248 may beone or more microprocessors, micro-controllers, field programmable gatearrays, or ASICs, while memory 250 may include random access memory(RAM), dynamic random access memory (DRAM), static random access memory(SRAM), flash memory, and/or another digital storage medium, and alsotypically implemented using circuit logic disposed on one or morephysical integrated circuit devices, or chips. As such, memory 250 maybe considered to include memory storage physically located elsewhere inthe general purpose processor 246, e.g., any cache memory in the atleast one CPU 248, as well as any storage capacity used as a virtualmemory, e.g., as stored on a mass storage device 252, or anothercomputing system through at least one network interface 256 (illustratedas, and hereinafter, “network I/F” 256) by way of at least one network258. It will be appreciated that the at least one network 258 mayinclude at least one private communications network (e.g., such as anintranet) and/or at least one public communications network (e.g., suchas the Internet). Similarly to the general purpose processor 246,computing system 254, in specific embodiments, is a computer, computersystem, computing device, server, disk array, or programmable devicesuch as a multi-user computer, a single-user computer, a handheldcomputing device, a networked device (including a computer in a clusterconfiguration), a mobile telecommunications device, a video game console(or other gaming system), etc.

The general purpose processor 246 is coupled to at least one peripheraldevice through an input/output device interface 260 (illustrated as, andhereinafter, “I/O I/F” 260). In particular, the general purposeprocessor 246 receives data from a user through at least one userinterface 262 (including, for example, a keyboard, mouse, a microphone,and/or other user interface) and/or outputs data to the user through atleast one output device 264 (including, for example, a display,speakers, a printer, and/or another output device). Moreover, in someembodiments, the I/O I/F 260 communicates with a device that isoperative as a user interface 262 and output device 264 in combination,such as a touch screen display (not shown). Furthermore, general purposeprocessor 246 may communicate through the Network IF 256 with RF frontend hardware 266, such as a downconverter and sampler, as illustrated.Alternately, the RF front end hardware 266 may communicate via the I/OI/F 260.

The general purpose processor 246 is typically under the control of anoperating system 268 and executes or otherwise relies upon variouscomputer software applications, sequences of operations, components,programs, files, objects, modules, etc., consistent with embodiments ofthe invention. In specific embodiments, the general purpose processor246 executes or otherwise relies on an application 270 to implement atleast a portion of the signal decomposition and parameterizationalgorithm to obtain estimates for direct and indirect rays present inreceive GNSS data consistent with embodiments of the invention.Moreover, and in specific embodiments, the general purpose processor 246may be configured with a ray database 272 to store estimate parametersfor use in inter-PIT algorithms consistent with embodiments of theinvention.

While the present invention has been illustrated by a description of oneor more embodiments thereof and while these embodiments have beendescribed in considerable detail, they are not intended to restrict orin any way limit the scope of the appended claims to such detail.Additional advantages and modifications will readily appear to thoseskilled in the art. The invention in its broader aspects is thereforenot limited to the specific details, representative apparatus andmethod, and illustrative examples shown and described. Accordingly,departures may be made from such details without departing from thescope of the general inventive concept.

What is claimed is:
 1. A method for intra-PIT signal decomposition of asignal received with RF front end hardware, the method comprising:aligning the signal into integer multiples of a duration of apseudorandom noise code sequence; computing a search grid based on thealigned signal; coarsely estimating a plurality of initial rayparameters associated with the computed search grid; using the coarselyestimated plurality of initial ray parameters to initiate fineestimation of the plurality of initial ray parameters utilizingstochastic search and optimization techniques; computing an estimateerror for the fine estimation of the plurality of initial rayparameters; and in response to determining a reduction in the computedestimate error, removing the fine estimation of the plurality of initialray parameters from the computed search grid.
 2. The method of claim 1,wherein the plurality of initial ray parameters is selected from a groupconsisting of: amplitude, carrier frequency, Doppler offset, propagationdelay, carrier phase, and combinations thereof.
 3. The method of claim1, wherein the signal received with the RF font end hardware is sampled,and wherein aligning the signal comprises: dividing the sampled signalinto vectors.
 4. The method of claim 3, wherein the sampled signal isdivided into PIT-sized vectors.
 5. The method of claim 3, whereincoarsely estimating a plurality of initial ray parameters comprises:estimating Doppler frequency offset, propagation delay, and carrierphase values; and calculating initial amplitude estimate using theestimated Doppler frequency offset, propagation delay, and carrier phasevalues.
 6. The method of claim 5, wherein initiating fine estimation ofthe plurality of initial ray parameters comprises: decomposing thecalculated search space; and computing an estimate search space.
 7. Themethod of claim 6, wherein decomposing the calculated search spacecomprises: computing search grid; constructing trial code delay vector;constructing trial ray waveform vector; computing trial estimate searchgrid using the trial code delay vector and trial ray waveform vector;and comparing received data search grid and computed trial estimatesearch grid.
 8. The method of claim 1, wherein the received signal is aGPS L1 C/A-coded signal.
 9. An apparatus, comprising: RF front endhardware configured to receive, downconvert, and sample a signal; amemory; a processor; and a program code resident in the memory andconfigured to be executed by the processor for intra-PIT signaldecomposition of the signal received by the RF front end hardware, theprogram code further configured to align the signal into integermultiples of a duration of a pseudorandom noise code sequence, compute asearch grid based on an integer multiple of the aligned signal, coarselyestimate a plurality of initial ray parameters associated with thecomputed search grid, use the coarsely estimated plurality of initialray parameters to initiate fine estimation of the plurality of initialray parameters utilizing stochastic search and optimization techniques,compute an estimate error for the fine estimation of the plurality ofinitial ray parameters, compute a stopping criteria statistic bycomparing a peak power of the search grid with a noise power present inthe search grid, and in response to determining the stopping criteriastatistic is less than a stopping criteria threshold, processing a nextinteger multiple of the aligned signal.
 10. The apparatus of claim 9,wherein the plurality of initial ray parameters is selected from a groupconsisting of: amplitude, carrier frequency, Doppler offset, propagationdelay, carrier phase, and combinations thereof.
 11. The apparatus ofclaim 9, wherein the program code is further configured to align thesignal by: dividing the sampled signal into vectors.
 12. The apparatusof claim 11, wherein the program code is further configured to dividethe sampled signal into PIT-sized vectors.
 13. The apparatus of claim11, wherein the program code is further configured to coarsely estimatea plurality of initial ray parameters by: estimating Doppler frequencyoffset, propagation delay, and carrier phase values; and calculatinginitial amplitude estimate using the estimated Doppler frequency offset,propagation delay, and carrier phase values.
 14. The apparatus of claim13, wherein the program code is further configured to initiate fineestimation of the plurality of initial ray parameters by: decomposingthe calculated search space; and computing an estimate search space. 15.The apparatus of claim 14, wherein the program code is furtherconfigured to decompose the calculated search space by: computing searchgrid; constructing trial code delay vector; constructing trial raywaveform vector; computing trial estimate search grid using the trialcode delay vector and trial ray waveform vector; and comparing receiveddata search grid and computed trial estimate search grid.
 16. Theapparatus of claim 9, wherein the signal received by the RF front endhardware is a GPS L1 C/A-coded signal.
 17. The apparatus of claim 9,wherein the stopping criteria statistic is computed using a cumulativedistribution function of a Rayleigh distribution.
 18. The apparatus ofclaim 17, wherein the cumulative distribution function is expressed as${F_{p_{ɛ}} = {1 - {\exp \left( {- \frac{\xi_{p_{ɛ}}^{2}}{2\sigma_{x}^{2}}} \right)}}},{\xi \geq 0.}$19. A method for intra-PIT signal decomposition of a signal receivedwith RF front end hardware, the method comprising: aligning the signalinto integer multiples of a duration of a pseudorandom noise codesequence; computing a search grid based on an integer multiple of thealigned signal; coarsely estimating a plurality of initial rayparameters associated with the computed search grid; using the coarselyestimated plurality of initial ray parameters to initiate fineestimation of the plurality of initial ray parameters utilizingstochastic search and simulated annealing; computing a stopping criteriastatistic by comparing a peak power of the search grid with a noisepower present in the search grid; and in response to determining thestopping criteria statistic is less than a stopping criteria threshold,processing a next integer multiple of the aligned signal.
 20. The methodof claim 19, wherein the plurality of initial ray parameters is selectedfrom a group consisting of: amplitude, carrier frequency, Doppleroffset, propagation delay, carrier phase, and combinations thereof. 21.The method of claim 19, wherein the signal received with the RF font endhardware is sampled, and wherein aligning the signal comprises: dividingthe sampled signal into vectors.
 22. The method of claim 21, wherein thesampled signal is divided into PIT-sized vectors.
 23. The method ofclaim 21, wherein coarsely estimating a plurality of initial rayparameters comprises: estimating Doppler frequency offset, propagationdelay, and carrier phase values; and calculating initial amplitudeestimate using the estimated Doppler frequency offset, propagationdelay, and carrier phase values.
 24. The method of claim 23, whereininitiating fine estimation of the plurality of initial ray parameterscomprises: decomposing the calculated search space; and computing anestimate search space.
 25. The method of claim 24, wherein decomposingthe calculated search space comprises: computing search grid;constructing trial code delay vector; constructing trial ray waveformvector; computing trial estimate search grid using the trial code delayvector and trial ray waveform vector; and comparing received data searchgrid and computed trial estimate search grid.
 26. The method of claim19, wherein the stopping criteria statistic is computed using acumulative distribution function of a Rayleigh distribution.
 27. Themethod of claim 26, wherein the cumulative distribution function isexpressed as${F_{p_{ɛ}} = {1 - {\exp \left( {- \frac{\xi_{p_{ɛ}}^{2}}{2\sigma_{x}^{2}}} \right)}}},{\xi \geq 0.}$